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ABSTRACT 


Flight experiments were defined for the COLD-SAT test bed satellite and the Shuttle 
middeck to help establish the influence of the gravitational environment on liquid slosh dynamics 
and control. Several analytical and experimental studies were also conducted to support the 
experiments and to help understand the anticipated results. 

Both FLOW-3D and NASA-VOF3D computer codes were utilized to simulate low Bond 
number, small amplitude sloshing, for which the motions are dominated by surface forces; it was 
found that neither code provided a satisfactory simulation. Thus, a new analysis of low Bond 
number sloshing was formulated, using an integral minimization technique that will allow the 
assumptions made about surface physics phenomena to be modified easily when better knowledge 
becomes available from flight experiments. Several examples were computed by the innovative 
use of a finite-element structural code. An existing spherical-pendulum analogy of nonlinear, rotary 
sloshing was also modified for easier use and extended to low-gravity conditions. 

Laboratory experiments were conducted to determine the requirements for liquid-vapor 
interface sensors as a method of resolving liquid surface motions in flight experiments. The 
feasibility of measuring the small slosh forces anticipated in flight experiments was also investigated. 
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1.0 INTRODUCTION 

1.1 Background 

Many of the future spacecraft planned by NASA and DOD will contain large quantities of 
liquids for propellants, life support, and other uses. As two examples, space vehicles for the 
ambitious Space Exploration Initiative will transport and store hundreds of thousands of kilograms 
of cryogenic propellants, and Space Station Freedom will contain tens of thousands of both cryogenic 
and storable liquids that must be re-supplied in orbit Future deep space probes and space-based 
optical systems will also have a large fraction of their mass in the form of liquids. Controlling, 
pointing, or docking of such spacecraft is critically dependent upon understanding and managing 
the motion (sloshing) of the large masses of liquids in their tanks. 

When liquid body forces are large, the "high gravity" motions of a liquid contained in a 
tank are well understood, and analytical, numerical, and scale-model test methods have been well 
established to treat them [1,2, 3]. But, when liquid body forces are small, the "low gravity" motions 
of liquid free surfaces are not nearly so well understood because they are dominated by surface 
forces that are completely masked on earth by the much larger body forces. Note that the true 
indicator of "low" gravity is the Bond number Bo, defined as Bo = g^R ^/ p where g# is the effective 
gravity or settling acceleration, R 0 is a representative dimension of the tank such as diameter, and 
P is the specific surface tension of the liquid; Bond numbers of ten or less are generally considered 
to represent "reduced" gravity and Bo « 1 to represent "micro" gravity. Spacecraft in near-earth 
orbit usually have Bo ~ 0.1 while deep space probes have Bo « 1. 

The rudimentary flight tests that have been conducted to date [e.g., 4] have provided little 
quantitative data about low gravity sloshing. Some information is available from drop tower and 
ground test simulations for Bo ~ 10 [e.g., 5, 6, 7] but these kinds of studies are hampered by the 
small size of the tanks that must be used, the short test time, and (at least up to the present) the need 
to use non-cryogenic liquids. As an example of the available information, it appears that viscous 
slosh damping £ is larger in low gravity and, for a cylindrical tank, can be correlated by [5, 7]: 

C = A{vlf K Rl) ia n+C{Bo)- b ] (1.1-1) 

where v is the kinematic viscosity, f n is the slosh natural frequency, and A, b, and C are 
empirically-determined constants. Although £ should increase in low gravity because of the larger 
wetted area of the liquid on the tank walls, the increase predicted by Eq. (1.1-1) for small Bo is 
much larger than can be accounted for merely by the larger wetted area. Since the low gravity 
response of a liquid to tank motions depends strongly on the physics of the contact line of the liquid 
at the tank wall [8,9], the anomalously high damping may be the result of poorly understood effects 
at the contact line. Furthermore, the natural frequency measured in ground simulations of low 
gravity sloshing varies by more than a factor of two between the extreme limits of a "free" contact 
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line and a "stuck" contact line [7]. For cryogens, the contact line motion may also be affected by 
evaporation and condensation effects. Contact line effects are totally negligible when liquid body 
forces are large. Nonlinear sloshing also appears to be more prominent in low gravity [10, 11]. 
Finally, because of the lack of low gravity slosh test data to provide insight, all methods to date for 
predicting liquid free surface motions in low gravity have employed potentially unrealistic 
assumptions about the contact line motion; none of the assumptions has as yet been adequately 
validated by comparison to test data. 

In order to gain the required physical understanding about liquid motions in low gravity, 
experiments are needed in an actual low gravity environment, using tanks of at least moderate size 
and long test durations. The study documented in this report defines such a set of space-based 
experiments and summarizes the supporting analytical models and ground-based testing. The study 
initially concentrated on experiments using liquid hydrogen that were to be conducted with 
NASA-LeRC’s COLD-SAT test-bed satellite, but it was later extended to include the definition of 
experiments using non-cryogenic liquids that can be conducted in the Space Shuttle. 

1.2 Objective 

The objective of the program was to conduct research to establish the influence of the 
gravitational environment on liquid slosh dynamics and control, including analytical and 
experimental studies. 

1.3 Scope of Work 

The project work was conducted in four technical tasks described below and one 
administrative task. Succeeding sections of this Final Report discuss the accomplishments of each 
technical task in detail. 

1.3.1 Task I — Technology Requirements 

For this task, future and planned NASA and DOD missions were reviewed to determine 
typical ranges of parameters that are important for liquid dynamics and control. The parameters 
included: tank shape, size, and fill level; internal tank structure and anti-slosh devices; spacecraft 
maneuvers; and liquid thermophysical properties. 

1.3.2 Task II — Definition of Flight Experiments 

The efforts of this task were devoted to defining the specifications, instrumentation, and 
data requirements for space experimentation on liquid free surface dynamics. It was composed of 
three subtasks: 

1. Review the preliminary COLD-SAT slosh dynamics experimental requirements 
document and prepare an updated version, based on the conclusions from Task I. 

2. Define an alternative flight experiment for the Space Shuttle, using a non-cryogenic 
liquid. 
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3. Provide technical oversight to NASA as requested on the detailed development of 
all COLD-SAT experiments. 

1.3.3 Task III — Analytical Model Development 

The efforts of this task were devoted to developing analytical methods and models that 
would allow liquid free surface experimental results to be analyzed and extended. It was composed 
of four subtasks: 

1. Investigate methods of improving the representation of surface physics effects in 
FLOW-3D (conducted by Flow Science, Inc.). 

2. Investigate NASA-VOF3D as a method of simulating liquid slosh in a low gravity 
environment. 

3. Develop a linearized slosh model of low gravity sloshing and implement it using 
available finite-element computational technology. 

4. Extend an existing spherical pendulum model of nonlinear rotary sloshing to low 
gravity conditions. 

Subtasks 3 and 4 constituted the bulk of the effort for Task III. 

1.3.4 Task IV — Ground Experimentation 

The efforts of this task were devoted to ground tests in support of instrumentation 
development for the space experiments. It was composed of two subtasks: 

1. Investigate liquid-vapor interface sensors as a method of tracking moving free 
surfaces in a low gravity environment. 

2. Determine the specifications for load cells and accelerometers that could be used to 
measure low gravity slosh forces and moments. 
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2.0 TECHNOLOGY REQUIREMENTS 

The definitions of low-gravity liquid-motion flight experiments were based on the ranges 
of parameters that are anticipated for future NASA and DOD missions. The parameters selected 
for this study were: (1) tank shape, size, and filling level; (2) internal tank structure and anti-slosh 
devices; (3) low gravity environment; (4) vehicle motions; and (5) liquid properties. 

2.1 Data Sources 

Over twenty-five technical reports discussing proposed NASA and DOD missions were 
reviewed to obtain numerical estimates of the selected parameters; references [12 - 17] are 
representative examples of the reviewed reports. In addition, cognizant personnel at NASA centers 
and at the Air Force Astronautics Laboratory were interviewed by project personnel. 

2.2 Representative Tank Shapes, Volumes, and Gravity Levels 

As specific examples, some details of the tanks for six planned or designed space vehicles 
are summarized below. 

Orbital Transfer Vehicle — The propellants are LH 2 and L0 2 . The baseline design LH 2 
tank is cylindrical with ellipsoidal ends and has a volume of approximately 5000 ft 3 ( 1 40 m 3 ) 
and a diameter of 14 ft (4.4 m). The L0 2 tank is a sphere with a volume of approximately 
1500 ft 3 (42 m 3 ) and a diameter of approximately 14 ft (4.4 m). The gravity environment 
ranges from about lg 0 to about 10' 5 g o . The tanks for the proposed European Space Tug 
are of similar sizes. 

CRAF/Casslnl Space Probe — The propellants for the main engine are monomethyl 
hydrazine and N 2 0 4 ; the propellant for the RCS engines is hydrazine. The main engine 
tanks are spheres with a diameter of 5.1 ft (1.56 m). The RCS tanks are cylinders with 
rounded ends, with a total length of 2.6 ft (0.8 m) and a diameter of 1.8 ft (0.56 m). The 
gravity environment ranges from about 0.01 g 0 during a main engine firing to practically 
zero gravity during a deep space coast. 

Liquid Droplet Radiator System — The planned tanks are of various shapes and sizes, 
depending on the design, but are typically spheres about 3 ft (1 m) in diameter. The tanks 
contain liquid metal or a very low vapor pressure oil. The gravity environment is lO^g,, 
to 10' 6 g o . 

Strategic Defense Initiative Optical Systems — Specific tank designs are classified. 
Typical tank shapes are spheres and cylinders with rounded ends. The liquids are cryogenic 
propellants as well as storables. One unclassified design study tank is a 21 ft (6.4 m) sphere 
containing up to 21000 lbs (9500 kg) of LH 2 . The gravity environment can be as low as 

io- 7 g 0 . 
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Compact Cryogenic Feed System — The LO z tank of this system is a torus with a 
major diameter of 9.6 ft (2.9 m) and a minor diameter of 4.6 ft (1.4 m). The gravity 
environment ranges from high-gravity to 10'^. 

Space Station — Some designs of the Space Station include tanks that are large enough 
to refill two Orbital Transfer Vehicle LH 2 tanks. Typical "large" tanks are spheres with 
volumes of 5000 ft 3 (140 m 3 ). The gravity environment is about 10‘ 5 g o - 10' 6 g o . 

2.3 Summary of Parameters 

Based on the data sources and interviews, it was concluded that the most common tank 
shapes for future missions are spheres and cylinders with ellipsoidal ends. Tank volumes range 
from 5 ft 3 (0. 14 m 3 ) to over 5000 ft 3 (140 m 3 ). During most missions, the tank filling level will vary 
from nearly full to nearly empty. Most tanks will contain liquid acquisition devices, typically 
screened channels running the length of the tank. None of the tanks surveyed contained any specific 
anti-slosh devices since the designs were not yet sufficiently mature to consider sloshing; however, 
liquid acquisition devices such as channels and central cruciform vanes are known to damp 
high-gravity sloshing effectively. 

The gravity environment covers the range from lg,, to lO^g,, (and even smaller for deep 
space probes). The non-dimensional Bond number of course depends on tank size and liquid surface 
tension as well as gravity level. As examples, at the lower range of gravity levels (*10‘ 6 g o ) the 
Bond number for LH 2 is 0.09 for a 1 m diameter tank and 9.1 for a 5 m diameter tank; for L0 2 , the 
corresponding Bond numbers are 0.21 and 21. The contact angle for most cryogens and propellants 
is close to 0* against tank materials of aerospace interest; that is, the liquid "wets" the tank wall 
material. The ratio P of surface tension to density for most cryogens and propellants ranges from 
0.0007 ft 3 /sec 2 (20 cm 3 /sec 2 ) to 0.001 ftVsec 2 (30 cm 3 /sec 2 ). 

Space vehicle motions are of three general kinds: (1) g-jitter; (2) attitude control maneuvers; 
and (3) thrusting and other large impulses such as docking. The magnitude and frequency content 
of g-jitter depends on the space vehicle in question; for the Space Shuttle, as an example, the rms 
magnitude is of the order of 10^g o and almost all the frequency content is above 1 Hz. Attitude 
control maneuvers are likewise dependent upon the mission; typical maneuvers include slewing 
around one of the vehicle axes at 2*/sec and low-frequency, small-amplitude oscillations around 
all three axes. Thrusting and docking can impose impulsive accelerations typically of the magnitude 
of 0.1 g c . The variety of tank motions implies that both large and small amplitude liquid motions 
will be excited. 

2.4 Conclusions 

Based upon the current state of knowledge and the survey described above, the desired 
technology requirements for a flight experiment can be summarized as: 


6 



• Spherical and cylindrical tanks with ellipsoidal ends should be tested; "bare" tanks 
and tanks with internal liquid acquisition devices are both important 

• Cryogenic liquids are desirable as the test liquid; in any case, the test liquid should 
wet the tank wall material. 

• Bond numbers for the tests should cover a range from =10 to near zero. 

• Tank motions should include large and small amplitude impulsive accelerations and 
sustained small-amplitude oscillatory accelerations. 

• Liquid motions investigated in the tests should include small amplitude, linear 
sloshing; moderate amplitude nonlinear planar and rotary sloshing; and large 
amplitude, reorientation-like motions that eventually settle down to long-lived 
small-amplitude sloshing. 

Data acquired from the tests should be sufficient in quality and quantity to validate analytical models 
and guide the development of improved models. Previous ground-based slosh tests have indicated 
that the most useful forms of data are (in descending order of usefulness): detailed maps of free 
surface shapes and motions (motion pictures, video recordings, etc.); slosh natural frequencies; time 
history or frequency sweeps of the forces and torques exerted on the tank; and slosh damping. 


7 



8 



3.0 DEFINITION OF FLIGHT EXPERIMENTS 

3.1 General 

A flight experiment was defined for both the COLD-SAT satellite and the Space Shuttle 
middeck. The complete Experiment Requirements Document for the COLD-SAT slosh dynamics 
experiment is included as Appendix A of this Final Report; it is summarized in Section 3.2. The 
definition of the Space Shuttle experiment is given in Section 3.3. The characteristics common to 
both definitions are described below. 

3.1 .1 Objectives 

The specific objective of both flight experiments is to quantify the liquid motion resulting 
from typical maneuvers of space vehicles in a low gravity environment. Physical processes to be 
investigated include: 

Static liquid configurations — The shape and location of the liquid in the test tanks 
will be monitored under ambient orbital conditions. 

Liquid response to various discrete accelerations — Impulsive and periodic 
accelerations of selected amplitude, frequency, and duration will be applied and the free 
surface response monitored. 

Slosh damping — Viscous damping will be determined for all tested conditions; if 
possible, the damping provided by a ring baffle should also be measured as a demonstration 
of an anti-slosh device for low gravity sloshing. 

3.1 .2 Key Parameters 

The key parameters that are expected to influence sloshing in low gravity are: (1) liquid 
surface tension, density, and viscosity; (2) contact angle of the liquid at the tank wall; (3) tank shape 
and internal hardware; (4) liquid fill level; (5) steady settling acceleration; and (6) disturbance 
acceleration. Because most liquids of aerospace technological importance have static contact angles 
of = 0* against common metals, the static contact angle parameter will not be varied in the flight 
tests. The steady settling acceleration can be varied in the COLD-SAT experiment alternative but 
not in the Shuttle experiment. All other parameters will be varied in the relevant non-dimensional 
form. 

3.1 .3 Measurements 

The free surface configuration, slosh force, and tank acceleration environment constitute 
the bulk of the test measurements. 

Free surface configuration — Mapping the static and dynamic configuration of the 
free surface is one of the most important methods of understanding and correlating 
experimental results. For example, the wave shape near the wall indicates whether the 
contact angle remains constant during sloshing ("free" contact line) or changes. Likewise, 
rotary slosh is most easily detected by observing the free surface motion. Ideally, the 
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mapping should employ both visual records such as motion pictures and an array of digital 
liquid-vapor sensors to track the position of the free surface quantitatively at a number of 
discrete locations. For the COLD-SAT experiment, visual observations are not possible, 
so liquid-vapor sensors must be used exclusively. 

Liquid-vapor sensors should be positioned in at least two perpendicular arrays, as 
illustrated in Figure 3.1-1. Ground testing (see Section 5.0) indicated that the overall free 
surface shape and motion can be resolved satisfactorily by a radial spacing of the sensors 
as large as 40% of the free surface radius and a vertical spacing as large as 10% of the free 
surface radius. The level of resolution provided by this spacing is satisfactory to determine 
the predictive accuracy of existing models. More fundamental information may be required 
to determine weak assumptions and potential sources of improvement in the models. For 
example, the relation between the dynamic contact angle and the contact line velocity can 
have a significant effect on the slosh frequency and force [7], even though existing models 
assume that the dynamic contact angle remains equal to the static contact angle. To acquire 
detailed data about the wave shape near the wall from which the behavior of the dynamic 
contact angle can be inferred, a denser array of sensors is required in the vicinity of the 
wall, for at least one liquid level. Sensors in these denser arrays would have to be spaced 
symmetrically at about 0.02R o above and below the selected free surface level over a total 
vertical distance of about 0.1R o , and at least two such vertically-spaced arrays would have 
to be positioned radially at intervals of 0.02R o from the wall. Because of the "bent over" 
geometry of the tank-liquid intersection for spherical tanks (see Section 4.0 and Figure 
4.4-3), the dense arrays are most readily made applicable to cylindrical (straight wall) tanks, 
although a dense array of sensors positioned along a radial line could be used for spherical 
tanks. In addition, the array installation must not significantly interfere with the slosh wave 
motion. It is understood that the use of dense arrays of liquid-vapor sensors will require 
substantial data acquisition rates and may interfere with other experiments; thus, their use 
may be restricted to the Shuttle-based experiment. 

Load cells — The slosh force and torque exerted on the tank by the liquid in motion will 
be measured by load cells. The support structure of the tanks must be designed to 
accommodate these cells. 

Acceleration environment — Both the effective settling acceleration and the time history 
of the disturbance accelerations will be measured by accelerometers and gyros. 
Temperature and pressure — Liquid temperature will be measured at the start and end 
of each test to determine the liquid properties. For the COLD-SAT experiment in which 
cryogens are used, liquid pressure will also be measured. 
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3.1.4 Data Analysis 

The initial configuration of the liquid and the static contact angle will be determined by 
liquid volume and liquid-sensor measurements (and visual records for the Shuttle experiment) before 
the perturbation acceleration is applied. Accelerometer and gyro data will be used to determine the 
motion imposed on the tanks. 

For the impulsive perturbation tests, liquid-vapor sensor measurements (and video or 
cinema recordings for the Shuttle experiment) will be analyzed to determine: (1) slosh wave shape 
and amplitude as a function of time; (2) slosh natural frequency; and (3) slosh damping from the 
decay of the slosh wave amplitude. If a dense array of sensors is used near the wall as suggested 
above, contact angle and contact line velocity will also be computed as a function of time; if not, 
these quantities will still be estimated but the resolution is not expected to be sufficient to identify 
the relation between dynamic contact angle and contact line velocity. Load cell force histories will 
be analyzed to confirm the slosh natural frequency and damping data and to compute the slosh 
force. For periodic perturbations, the data analysis will be similar to the impulsive acceleration 
tests. In addition, load cell data, in conjunction with damping data from the impulsive tests for the 
same Bond number and liquid level, will be analyzed to determine the effective mass of liquid 
participating in the sloshing. When nonlinear effects are prominent, the line of action and the 
phasing of the slosh force relative to the excitation will also be determined from the load cell data. 

Physical properties of the liquids will be computed from the temperature and pressure 
measurements, using tabulated data. 

3.2 Summary of COLD-SAT Experiment Definition 

3.2.1 General Description 

Baseline designs of the COLD-SAT satellite envision three LH 2 tanks; a large "supply" 
tank [volume = 125 - 175 ft 3 (3.5 - 5 m 3 )] and two "receiver" tanks, one that is almost spherical 
[diameter =3.0 ft (0.9 m)] and one that is cylindrical with ellipsoidal ends [diameter = 3.0 ft (0.9 m), 
length = 5.5 ft (1.7 m)]. Slosh tests are defined for both receiver tanks to investigate the influence 
of tank shape. The cylindrical receiver tank contains a liquid acquisition device (screened channels), 
and various spray systems to study tank filling, cooldown, and cryogenic mixing. The spherical 
receiver tank is nominally "bare" but does contain spray systems. For the slosh tests, it is proposed 
that the cylindrical tank be fitted with a single ring baffle near its midpoint, to permit the investigation 
of a typical anti-slosh device. The COLD-SAT propulsion system will be used to provide the steady 
settling accelerations and the disturbance accelerations for the slosh tests. The slosh experiments 
can be performed as opportunities permit when the liquid filling levels of interest are available 
during other tests. 
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3.2.2 Measurements and Instrumentation 

The measurement and instrumentation requirements are summarized in Table 3.2-1. 

3.2.3 Hardware Requirements 

The special hardware items required for the experiment are: (1) annular ring mounted in 
the cylindrical receiver tank; (2) load cells and accelerometers mounted on the support structures 
of both receiver tanks; and (3) liquid-vapor sensor arrays in both receiver tanks. 

3.2.4 Test Matrix 

Two general kinds of tests will be conducted: (1) oscillatory disturbance acceleration, and 
(2) impulsive disturbance acceleration. The tests are further arranged as being of first and second 
priority. Table 3.2-2 summarizes the test matrix. The steady settling accelerations listed in the 
table conform to the nominal levels that can be obtained by firing a single COLD-SAT engine or 
a combination of engines. The nominal perturbation accelerations correspond to levels that can be 
obtained by firing the the appropriate attitude control thrusters or, depending on the COLD-SAT 
design, a gimbaled engine; impulsive accelerations are obtained by firing the engine for a short 
time and periodic accelerations by on-off firings for scheduled periods. The predicted natural period 
of the slosh waves for the various tests is indicated by the symbol x, under the heading "Comments" 
in the table. 

3.3 Shuttle Flight Experiment 

3.3.1 General Description 

The Shuttle flight experiment is designed to fit into two middeck lockers using a double 
mounting plate. Figure 3.3-1 shows a conceptual design of the flight experiment package. Three 
spherical tanks will be used in the first flight; cylindrical tanks will be used for a second flight unless 
the tanks can be changed out during the first flight. Two alternatives are considered: all tanks have 
different diameters [13 in., 7 in., and 3.5 in (33 cm, 18 cm, 9 cm)] to permit a thorough investigation 
of tank size effects, or the two smaller tanks have the same size [7 in (18 cm)] to permit liquids of 
different viscosity to be tested in tanks of the same size and shape. Further trade studies are required 
to determine the relative advantages of each alternative. The test liquid is either water (containing 
a surfactant, defoamer, and bactericide) or silicone oil, the choice depending primarily on safety 
trade studies. 

Figures 3.3-2a and 3.3-2b show the Bond number and predicted slosh frequency for a 
spherical tank as a function of settling acceleration and tank diameter, for a liquid having 
o/p = 0.0009 ft 3 /sec 2 (25 cm 3 /sec 2 ) typical of silicone oils. (The Bond numbers and frequencies 
would be slightly larger for water.) For all the tank sizes suggested for the flight experiment, the 
Bond number is considerably less than one for the steady settling accelerations that are anticipated 
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TABLE 3.2-1 MEASUREMENT REQUIREMENTS FOR COLD-SAT LIQUID DYNAMICS EXPERIMENT 
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TABLE 3.2-2 COLD-SAT TEST MATRIX 
First Priority Tests 

Tank Fill Axial g Bond Perturbation Test Time Comments 
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for the Shuttle. The flight experiment does not contemplate thrusting to increase the settling 
acceleration to give larger Bond numbers. The use of very small Bond numbers has the 
advantage that the slosh characteristics for a given tank size and shape are nearly independent of 
the magnitude of the setding acceleration; hence, accurate measurements of the Shuttle ambient 
acceleration are not required. All the tests will be conducted with half-full tanks that have been 
pre-loaded and sealed on the ground. 

The desired perturbation accelerations will be provided by a controllable electrodynamic 
shaker attached between the locker mounting plate and tank support plate, as shown in Figure 3.3-1. 
When the test tanks are cylindrical, the axis of the tanks must be aligned with the steady g-vector 
at the locker location to within ± 2* to preserve the symmetry of the liquid free surface configuration. 
To accomplish this while allowing some flexibility in the experiment integration and mission 
planning, the orientation of the shaker and tank support plate combination can be manually adjusted 
over a range of about ±15* before the experiment is installed in the locker. Orientation of the 
spherical test tanks with respect to the g-vector is not critical, but the axis of the electrodynamic 
shaker should still be roughly perpendicular to the g-vector in order to excite lateral sloshing. 

For the proposed tank sizes, the liquid motions of the types to be studied in the experiment 
are insensitive to vernier RCS firings because the magnitude and duration of the resulting tank 
accelerations are less than that required to destabilize the free surface. They should also be 
insensitive to g-jitter because of the large mismatch between the g-jitter frequencies and the slosh 
natural frequencies. Accelerations from a main RCS firing will, however, cause the liquid surface 
to break up or grossly reorient, so the longer duration tests must be conducted during quiet times. 

The weights, power, and size requirements are within middeck allowables, assuming that 
a Shuttle video camcorder can be supplied by NASA without penalty to the experiment. Trade 
studies during further analysis (i.e., Phase B of a normal flight experiment program) may indicate 
that motion pictures offer significant advantages over video recordings; if so, the experiment can 
be expanded to three middeck lockers to permit additional battery power for the motion picture 
cameras and lighting. The c.g. location of the flight experiment package is about 10 in. (25 cm) 
from the face of the double adapter plate. Lockout mechanisms will be installed between the tanks 
and the support plate and between the support plate and the locker side braces to eliminate potential 
damage during launch and re-entry. 

3.3.2 Measurements and Instrumentation 

The static orientation and the dynamic motions of the liquid free surface will be measured 
by liquid-vapor sensors and recorded by a camcorder or motion picture camera focused on the test 
tank. Slosh force and torques will be measured by load cells. The tank table accelerations will be 
measured by accelerometers mounted on the table. The table displacement and frequency will be 
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measured by a linear variable differential transformer (LVDT). All the measurements are within 
the capabilities of instrumentation used in previous ground studies. Figure 3.3-3 shows a block 
diagram of the instrumentation. 

A controller for the shaker actuator is included in the design to produce the desired excitation 
waveforms. Frequencies of 0.003 to 0.1 Hz are required. An onboard data logger - controller - 
mass storage computer will be used to acquire up to 60 channels of analog data which will be 
converted to 10-bit digital form for storage. All the instrumentation and data acquisition equipment 
will operate on the 28 volt d.c. power available at the middeck lockers. 

3.3.3 Hardware Requirements 

Table 3.3-1 lists the weights, sizes, and power requirements of the primary subsystems of 
the flight experiment. The listed weights, power consumptions, etc., are based on commercially 
available equipment. 

3.3.4 Test Matrix 

The test matrix is shown in Table 3.3-2. Five types of tests are defined: 

Ambient configuration — The liquid interface shape will be measured before each test 
(i.e., without excitation). The measurements will be used to validate analytical models and 
to confirm the desired initial conditions. 

Small amplitude Impulsive translation — An impulsive displacement will be applied 
by the shaker. The free decay of the wave motion (force and wave amplitude) will be 
recorded for up to eight slosh cycles. The tests will establish the fundamental slosh 
frequency, the frequency of any excited higher modes, and the viscous damping of the 
fundamental mode, all of which can be used to validate analytical models. Additionally, 
the results will allow the controller software to be updated to make subsequent harmonic 
excitation tests less time consuming. 

Small amplitude harmonic translation — The table frequency will be swept from 25% 
below the fundamental slosh frequency to 25% above it. Force amplitude and phase angle 
and surface wave amplitude will be measured continuously. These tests can be used to 
improve and validate analytical models and to establish equivalent mechanical models of 
low gravity sloshing. 

Large amplitude harmonic translation — The table frequency will be held constant 
for 8 to 10 cycles at a few specific frequencies near the slosh natural frequency. The 
amplitude of the shaker actuator will be considerably larger than for the previously 
described tests in order to excite nonlinear, rotary sloshing. Slosh force (in-line and 
cross-axis) and phase angle and free surface orientation will be measured. The results will 
be used to assess the importance of nonlinear effects in low gravity and to validate and 
improve analytical models. 
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TABLE 3.3-1 . SUBSYSTEM POWER REQUIREMENTS, SIZES, AND WEIGHTS 

FOR SHUTTLE FLIGHT EXPERIMENT 
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Large amplitude Impulsive translation — A large translational impulse will be 
applied to all the tanks simultaneously. After the initial large liquid motion is over, the 
setding time will be determined by measuring free surface orientations. The results can 
be used to improve and validate computational fluid dynamics codes of liquid 
reorientation and setding. 

The total test time, not necessarily continuous, is about four hours. Individual tests, which must be 
continuous, last from about 7 minutes to 35 minutes. 
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4.0 ANALYTICAL MODELS 

4.1 General 

Several different studies were pursued for prediction of liquid motions in low gravity. First, 
two computational fluid dynamics codes that can model free surface motions — FLOW-3D and 
NASA-VOF3D — were investigated to determine their ability to simulate sloshing at small Bond 
numbers. Both codes are widely used to simulate liquid motions in low gravity when inertia forces 
control the motions, as, for example, occurs when a large acceleration is induced by thrusting. The 
present investigation indicated, however, that neither code yields a satisfactory simulation of liquid 
motions when surface forces are dominant, such as might occur during spacecraft station-keeping 
or guidance maneuvers or after large amplitude motions have decayed to small amplitude sloshing. 
For that reason, an analysis of small-amplitude low gravity sloshing in axisymmetric tanks was 
developed here and implemented computationally using available finite-element technology. This 
model, which is described in Section 4.4, is meant to supplement FLOW-3D or NASA-VOF3D. 
In addition, an existing spherical pendulum mechanical analogy of nonlinear, rotary sloshing was 
modified to allow greater ease of numerical predictions and extended to low gravity conditions. 
This development is described in Section 4.5. 

4.2 FLOW-3 D Simulations 

FLOW-3D (1988 version) was used to simulate sloshing in a one-meter spherical tank for 
Bond numbers of 0.1, 1.0, and 10.0. The contact angle for the simulations was 5* and it was held 
constant during motion by activating the "wall adhesion" option. The equilibrium shape and position 
of the free surface were computed separately and input to FLOW-3D as an initial condition. A 
small lateral velocity was imparted to the tank to initiate the sloshing. 

The simulation was poor for all Bond numbers. In fact, the gas bubble above the free 
surface did not remain at its equilibrium position even when no motion was imparted to the tank, 
but instead migrated to the center of the tank. This behavior prevented the computation of realistic 
sloshing motions. After discussing the results with Flow Science, Inc. (the developers of 
FLOW-3D), it was concluded that the way surface tension and wall adhesion are modeled in 
FLOW-3D is not adequate to represent conditions in which surface tension forces are large compared 
to body forces. 

Consequently, Flow Science, Inc., through a subcontract with SwRI, reviewed the modeling 
of surface effects in FLOW-3D and suggested several possible improvements that would in principle 
allow the code to be extended to small Bond number conditions; the suggested improvements are 
documented in [18]. As a result of this review, a significant improvement in the wall adhesion 
model was incorporated in the 1990 version of FLOW-3D, which removes the tendency of the gas 
bubble to migrate away from the walls toward the center of the tank. A more exact representation 


21 


PRECEDING PAGE BLANK NOT FILMED 



of surface tension forces is also being evaluated but not yet incorporated in FLOW-3D. Until these 
improvements are incorporated, FLOW-3D is not capable of realistic simulations of sloshing when 
surface tension forces are predominant. 

4.3 NASA-VOF3D Simulations 

NASA-VOF3D and FLOW-3D both share the same structure and basic VOF (volume of 
fluid) algorithms. In particular, the representation of surface tension forces is similar for both codes, 
although NASA-VOF3D does not specifically contain a "wall adhesion" option to maintain a 
constant contact angle. Los Alamos National Laboratory (LANL) has modified the original VOF 
series of codes, under NASA sponsorship, to make them more applicable to fluid motions in low 
gravity. Since NASA-VOF3D is a vectorized code that executed quite slowly on SwRI’s VAX 
8700, SwRI requested LANL personnel to run the same three simulations attempted with FLOW-3D 
on LANL’s CRAY computer. The sequence of visualizations shown in Figure 4.3-1 gives some 
typical results from the LANL simulations. The bulk of the liquid surface does move in a 
sloshing-like wave motion but the contact line appears to be stuck to the wall. Furthermore, 
numerical convergence problems (which LANL has since eliminated) caused a fictitious wave to 
appear on the surface at r = 0 which then propagated away from the center. It was concluded that 
the surface force representation in NASA-VOF3D must be improved before realistic simulations 
of small Bond number sloshing are possible. 

Both FLOW-3D and NAS A-VOF3D solve unsteady flow problems by time-stepping from 
a known initial condition to a final condition that is compatible with a specified motion of the tank. 
They are thus inherently clumsy to use as a method of predicting, for example, the natural frequency 
of an oscillatory liquid motion. Consequently, as a tool for control system simulations or as a means 
of establishing equivalent mechanical models, neither code is as useful as the "eigenvalue" type of 
codes routinely used for high gravity sloshing. 

4.4 Linearized Low Gravity Sloshing Model 
4.4.1 Background 

Since the surface phenomena which can dominate low gravity motions are not yet well 
understood, previous analyses of low gravity liquid sloshing [7, 19-23, including FLOW-3D and 
NASA-VOF3D] have assumed simplified and perhaps unrealistic surface conditions; a typical 
assumption is that the contact angle remains constant during liquid motions, independent of the 
contact line velocity (i.e., the contact line is "free"). In the analysis described here, the equations 
of linear low-gravity sloshing are formulated and reduced to an integral-minimization technique 
that, while at present employing the same simplified surface physics assumptions as previous 
analyses, does permit a better representation of surface physics phenomena to be incorporated 
readily when the needed understanding becomes available. Errors in previous analyses of the low-g 
slosh force are also corrected in the present treatment. The equations are solved by a numerical 
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method based on a finite element structural code; the solution method can thus be easily adopted 
by other investigators. Several numerical examples are presented, and the parameters of an 
equivalent mechanical model are also computed. 

4.4.2 Equations of Motion 

Most tanks used in space applications are axisymmetric, and most liquids used in space 
applications are incompressible and nearly inviscid. Thus, the following idealizations are employed 
in the analysis: 

• the tank and static liquid configuration are axisymmetric; 

• the liquid motion is irrotational and derivable from a velocity potential <J>; 

• the motions are small enough to permit the equations to be linearized; and 

• only those slosh modes are considered that vary as the cosine of angular coordinate 0, 
since such modes are the only ones that produce a net force or moment on the tank. 

When needed, viscous damping can be incorporated as a small correction by the methods described 
by Abramson, et Al. [1]. 

The geometry is shown in Figure 4.4- 1 . Because the contact angle y c is small, the equilibrium 
free surface is highly curved. (Both the bottom and top of the tank can be dry for some combinations 
of low filling levels, tank shape, and y c [24], but such cases are not considered here.) A 
surface-normal coordinate system s,9,n is used to analyze the boundary conditions at the free 
surface in order to avoid the possibility of double valued expressions. 

The velocity potential must satisfy the condition of liquid incompressibility: 


v *~7Tr 


8<|> 

a 7 


1 3*$ a 2 * 

r 2 ae 2 dz 2 


=o 


in the fluid volume, V 


(4.4-1) 


Since the liquid cannot penetrate the tank walls, 4> must satisfy the "no flow" condition at the walls: 
_ _ a<t) 

V(b • n = — -0 on the walls, w (4.4-2) 

an 


where n is the outward-pointing normal. At the free surface, the wave velocity must be compatible 
with the liquid velocity: 

an ad) 

— = — on the free surface, / (4.4-3) 

at an 


where T) is the wave height measured normal to the free surface, as shown in Figure 4.4-1. 

The wave height and liquid velocity must be interrelated so as to satisfy the requirement 
of constant pressure p at the free surface. This relationship can be expressed as: 
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on the free surface, / 


(4.4-4) 
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where q(t ) is a function of at most time, and r(s) is evaluated on/. Since the liquid pressure at the 
surface is less than the gas pressure immediately above the liquid because of surface tension,/? must 
also satisfy the relation that: 

p =p g - 2gk on the free surface, / (4.4-5) 


where 2k is the mean curvature of the surface. From differential geometry, the mean curvature is 
related to the equation of the moving free surface, say ?(s,0,r) = n -ri(j,0,f) = O, by 
2k = -V • (V^ | | _1 ). Linearizing this expression for 2k with respect to the slosh wave height 

gives: 
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Here r h r 2 are the principal radii of curvature of the equilibrium free surface: 
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(4.4-6) 


(4.4-7) 


where r and z are evaluated on /. By combining Eqs. (4.4-4)-(4.4-7), subtracting out the 
time-independent terms (which are the equations of the equilibrium free surface described later), 
and absorbing q(t) into the gas pressure, the "dynamic" boundary condition of constant pressure at 
the free surface is obtained: 
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(4.4-8a) 


Equations (4.4-3) and (4.4-8a) can be combined to eliminate the wave height Tj: 
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<9r 2 

A free contact line (i.e., constant contact angle) is sometimes a reasonable approximation 
for some liquids and tank wall materials [7]. Hence, pending the acquisition of better knowledge 
from space experimentation, that condition is assumed here. With reference to the sketch shown 
in Figure 4.4-2a, this assumption reduces to: 


cosy, = n w (t) • «„(0 = (n w + An w ) • (n^ + An n ) = n w • n n + n w • A^ n + n n • An w (4.4-9) 


For the equilibrium surface, the unit vectors in Eq. (4.4-9) can be expressed as: 

n w = cos y c 7 H - sin y c 7 t ; «„ = ?„ (4.4- 10a) 

and to the first order in the arc lengths As and A£, the changes in the unit vectors as the surface 
moves away from its equilibrium position are: 
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(4.4- 10b) 


J(cosYA + sinY c eJ (4.4- 10c) 

Since Ay = rj cotY c and A£ = r|/ sin y c , the assumed condition at the contact line thus reduces to: 

= 0 at the contact line, s = s c (4.4- 11) 
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Note that Eq. (4.4-1 1) is satisfied identically when y e = 0*. Methods of incorporating contact line 
conditions other than the free condition will be discussed later. 


4.4.3 Equilibrium Surface 

Dropping the time-dependent terms from the combination of Eqs. (4.4-4) and (4.4-5) gives 
the equation of the equilibrium free surface z = f(s ): 


( 


1 1 




P * /- °U + ^ =Po ~ Pg=K 


(4.4-12) 


with boundary conditions: 

df 

0 at j = 0 (r = 0, z = —z 0 ) ; X“V = Yc at s =s c (r = r c , z =z e ) (4.4-13) 


The second part of Eq. (4.4-13) merely states that the slope angle of the free surface and the slope 
angle of the wall differ by the contact angle at the contact point. The pressure-reduction parameter 
Xo and the coordinates z 0 , s c , r c , and z c must be determined as part of the solution such that the 
volume of liquid under the free surface is equal to the specified liquid volume. 

4.4.4 Slosh Force 

In normal gravity, the slosh force is due entirely to the time-dependent liquid pressure 
exerted on the wetted walls. In low gravity, however, the "pull" of surface tension on the tank walls 
becomes unbalanced, which creates an additional force, and the pressure force itself includes a 
component that is negligible in normal gravity. 

The part of the lateral slosh force that is caused by the pressure on the tank walls is: 

/•2k /•/ c +V“X 

F p — [p 0 - P£z - p(9<]>/5r)]| r w cosQdzdB (4.4-14) 

In this equation, Tj, is the wave height at the contact line measured tangential to the wall (i.e., rj, 
corresponds to the wetted wall), and rj, sin y is the vertical displacement of the wave at the contact 
line (see Figure 4.4-2a). Note that T), = Ti c cotY;, where r) c =q (5 = s c , 0, r) is the wave height measured 
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normal to the free surface at the contact line. Since the slosh modes of interest vary with cos 0, the 
potential can be expressed as <J> = <])(r , z , t) cos 0 and the wave height as T] = T](,s,r) cos 0. Thus, Eq. 
(4.4-14) can be integrated with respect to 0 and then linearized with respect to <J) and r\ to give: 

F P = l(Po ~ Pgfc)^c cot Y c sin xln £ - Jtp £ 0^01 (4.4-15) 

The term p„ - p gf c , which is equal to -2 ok c , represents the pressure discontinuity from the gas to 
the liquid at the contact line that is caused by capillarity. Some previous analyses [e.g., 21] have 
linearized Eq. (4.4-14) before performing the integration and thus have incorrectly neglected this 
term. 

For the equilibrium surface, surface tension pulls on the tank walls equally around the 
circumference of the contact line, and there is clearly no net lateral force. But during sloshing, the 
contact line is displaced along the walls, and the surface tension pull becomes unbalanced. With 
reference to Figure 4.4-2b, a surface tension force cda is exerted on a differential element da of 
the displaced contact line. The component of this force in the plane of the wall is cdacosy e . The 
force is inclined to the horizontal at the small angle tan' , (-9r| f / r w dQ) ~ -5r\, / r w 90. Since da ~ 
r w d0 to first order, the net lateral force F a exerted on the tank by surface tension at the contact line 
is therefore: 

f 2 * — 

F a = -acosy.coty. (9r|/90)| r _ r sin 090 = (no cosy, cot yjTi,. (4.4-16) 

Jo ' 

With the exception of [7], previous analyses of the low-g slosh force have not considered F a . 

The total lateral slosh force F p + F a is: 

M — — 

F = — rcp J (d$/dt)r w dz - [7tocoty.(2K c r <: siny-cosy.)]^ (4.4-17) 

When the first part of Eq. (4.4-15) is neglected [21], the calculated slosh force will be too large. 
When F a is neglected [19, 20], the slosh force will be too small and can even be of the wrong sign. 

The slosh moment can be developed similarly but, for brevity, is not presented here. It is 
included, however, in the computer programs used to make numerical predictions. 

4.4.5 Non-Dimensional Equations 

For numerical work, it is best to solve the equations in non-dimensional form. Furthermore, 
since free vibrations are of interest, the equations are reduced to an eigenvalue problem by assuming 
that 4> varies harmonically in time as exp(icor), where to is the slosh natural frequency to be 
determined. Upper case letters are used for non-dimensional coordinates (for example, S = s/R 0 ). 
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The non-dimensional potential <I> is defined as 4>/V(l + Bo)oR 0 /p, the non-dimensional frequency 

Q as (0/^(1 + Bo)a/pRo, and the non-dimensional wave height H as iJ]/R 0 . After cancelling the 
common factor exp(t'cor) from all the equations, the non-dimensional form of the equations are: 
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(4.4-18) 

(4.4-19) 
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(4.4-22) 


The non-dimensional slosh force amplitude is: 
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cR 0 
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(4.4-23) 


where as before a bar over a quantity means that the 0 dependency has already been removed. 
Equation (4.4-23) clearly shows that the force components caused by surface tension (i.e., 
proportional to the wave height H c ) can be significant when the Bond number is small. 

4.4.6 Equivalent Mechanical Model 

The dynamics of linear sloshing can also be represented exactly by an equivalent mechanical 
model [1]. An appropriate model is a pendulum of length l, that has a mass m s which represents 
the liquid fraction that participates in the fundamental mode of the sloshing, and a rigidly-attached 
mass which represents the rest of the liquid. (The mass participating in the higher frequency sloshing 
modes is usually negligible.) However, in contrast to normal sloshing, the pendulum of the 
low-gravity model must be attached to the tank through a torsional spring k, which represents the 
stiffening effect of surface tension. From the non-dimensional process described in the previous 
section, the slosh natural frequency is expressed as: 

to 2 = Q 2 [(l+£o)a/ptf 0 3 ] (4.4-24) 


whereas the pendulum frequency is: 
<J pt nd = kAm,Z) + g/l, 


(4.4-25) 


Since these two frequencies must be the same for dynamic similarity, the parameters of the pendulum 
can be seen to be: 
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/. = RM 2 


(4.4-26a) 


K = Ctamjt/pRi = cR 2 m,/pR 3 0 Q 2 


(4.4-26b) 


Note that the ratio of the spring moment to the gravity moment, k,/mj,g, exerted on the pendulum 
pivot is equal to 1/Bo, and thus becomes negligible when Bo is large; that is, the mechanical model 
reduces to the usual high gravity form when Bo is large. 

The slosh mass m t can be computed from the amplitudes of the slosh force and liquid kinetic 
energy, as can be seen by considering the equivalent development for the mechanical model. For 
the model, the amplitude of the force when the pendulum is undergoing free oscillations is 
Fmod'i = m ,Q? x o> where x 0 is the amplitude of the oscillation, and twice the kinetic energy is 
2KEmod*i — m,(£> 2 x 2 . The slosh mass is therefore equal to the ratio of the force squared to twice the 
kinetic energy, a result that is independent of the arbitrary amplitude x a . The slosh mass can thus 
be computed by forming this ratio for the liquid. Twice the kinetic energy of the liquid is: 

2 KE = pJJ<j>(3<}>/3rt)<£4 (4.4-27) 

where the integral is evaluated on the free surface. The combination of Eq. (4.4-23) and the 
non-dimensional form of Eq. (4.4-27) then gives an expression for the slosh mass of the mechanical 
model: 
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(4.4-28) 


Although not developed here, an expression for the location of the pendulum pivot can be determined 
similarly from the slosh moment exerted on the tank. It is included, however, in the computer 
programs used to make numerical predictions. 

4.4.7 Integral Formulation 

Equations (4.4- 1 8)-(4.4-22) completely describe linear sloshing for the stated assumptions. 
They are, however, not solvable in closed form except when y c = 90* and then only for tank shapes 
that are coordinate surfaces. Hence, an approximate analytical method or a numerical method is 
needed. Most such methods apply only to the particular set of equations for which they have been 
developed and cannot be generalized. The solution method can be made much more amenable to 
changes in the surface physics assumptions by transforming the basic equations to an integral 
function form. The functions that minimize the integral are the desired solution. The appropriate 
integral can be derived by considering the kinetic, gravitational, and surface energies of the sloshing 
[25], or by partial integration and combination of Eqs. (4.4-18)-(4.4-22). The integral is: 
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/(<D,H) = i(l+So)aJJJv<J>* V<DrfV-(l+Bo)Q 2 JJ(<DH)dA 
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dAr 


(4.4-29) 
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This formulation has the additional advantage that, if trial solutions for H and <J> can be derived 
from any convenient method, the actual sloshing potential can be determined by using the trial 
functions to minimize /(<!>, H). The trial functions do not necessarily have to satisfy all the boundary 
conditions. The biggest advantage is, however, that an improved contact line condition can be 
treated by modifying the line integral in Eq. (4.4-29), without changing the solution method. 

4.4.8 Structural Finite-Element Simulation 

A finite element structural code was selected as the general method of finding trial solutions 
of Eqs. (4.4-18)-(4.4-22), for two reasons: (1) such codes can easily model the irregular and 
three-dimensional shapes that are typical of low-g liquid geometries, and (2) expertise in the use 
of such codes is widespread. 

By suppressing all the displacements of an elastic body but the x-component U, a structural 
finite element code can be forced to solve the following field equation [26]: 


( X + 2G ) 

d?U tfU d?U X+G 

{ tfu ?u) 

i 

{ G J 

dx 2 dy 2 dz 2 G 

^dxdy dxdz j 


G J * G dt 2 


(4.4-30) 


where X and G are elastic constants, J x is the body force, and p, is the material density. By arbitrarily 
choosing X = —G and lF x = p, = 0, Eq. (4.4-30) can be further reduced to: 


d 2 U tfU tfu 
■+— — -+- 


3 V 2 U = 0 


(4.4-31) 


dx 2 dy 2 d z 2 

which is the desired governing equation for potential flow, Eq. (4.4-1) or Eq. (4.4-18). 

A reactive traction force can be applied at any surface node (i.e., the nodes at the free 
surface or at the tank walls). The general form of this force is: 

Fj_ _ dU 

dn 


GA, 


dU 

a,U + a ->— — h 
dt 


tfU 1 


(4.4-32) 


where A f , is the surface area of the finite element in question. Equation (4.4-32) will be used to 
simulate the slosh boundary conditions. For nodes on the tank walls, the no- flow boundary condition 
can be treated exacdy by choosing a, = a 2 = a 3 = a 4 = 0 since this makes dU/dn = 0 at the nodes. 
However, the boundary condition of constant pressure at the free surface cannot be treated so easily. 
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The best feasible approximation to the free surface condition is to set a l = a 2 = a 4 = 0 and to choose 
a 3 (which is a mass-like parameter) as a representative value of all the terms within the square 
brackets in Eq. (4.4-8b); this choice makes dU/dn proportional to -dPU/dt 2 , which is at least of the 
same form as Eq. (4.4-8b). [Since the a, must be calculated in advance, they cannot depend on U, 
but this is what would be required to make Eq. (4.4-32) exactly the same as Eq. (4.4-8b).] The 
boundary condition that is actually applied to the nodes on the free surface is therefore: 


§ ■ 

or, after eliminating the exp(t'Qx) factor: 
cilJ 

— = (1 +Bo)C?M(S)U 
6N 


(4.4-33a) 


(4.4-33b) 


where a 3 = (1+Bo) M(S) is the specified distribution of fictitious masses at the surface nodes and 
the factor 1+Bo is included for analytical convenience. 

Using the values for the a, determined in this way, the structural code will yield a set of 
trial potential functions, say (/, for i = 1,2,3,...,K, and corresponding trial eigenvalues Cl 2 . The 
velocity potential for low-g sloshing is then expressed as: 

<D=I &,£/, (4.4-34) 

i = i 


where the b t are arbitrary parameters that must be chosen to minimize the integral, Eq. (4.4-29). 
Because the U i trial functions identically satisfy Eqs. (4.4-18)-(4.4-20), the minimization process 
is simplified considerably. In addition, M(S), as discussed below, will be chosen to satisfy the 
contact line condition, Eq. (4.4-1 1), which further simplifies the minimization process. The result 
is that Eq. (4.4-29) reduces to the simpler form: 



Bf ^ 

R dS 


as 
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+ ir R 


^ + y]w-^ i+so)f/ >)§ <is=o <4 - 4 - 35> 


for i = 1,2,3..., AT, evaluated on the free surface.. Integration with respect to 0 has already been 
performed. From Eq. (4.4-33), the dUfoN terms in Eq. (4.4-35) can be replaced by the corresponding 
(1 +Bo)MC%U , terms. The integration can therefore be performed readily using only the nodal 
values Ui along the line 0 = 0° (i.e., without any need for numerical differentiation of U t ). Equation 
(4.4-35) is a matrix eigenvalue problem which yields the non-dimensional slosh frequency £2 2 and 
eigenvectors ft,-. For cases when the trial functions cannot be made to satisfy the contact line 
condition, the contribution from the appropriate form of the line integral in Eq. (4.4-29) must also 
be included in Eq. (4.4-35). 
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4.4.9 Mass Distribution, M(S) 

In common with most integral-approximation solution techniques, it is relatively easy to 
insure that the numerical method predicts natural frequencies, or eigenvalues, with good accuracy. 
It is more difficult to insure that mode shapes, or eigenvectors, are predicted accurately. Nonetheless, 
it is important to predict the mode shapes accurately, because, as shown by Eq. (4.4-23), the slosh 
force, and therefore the slosh mass, is a function of the point value of the wave amplitude at the 
wall rather than on an integral average. Equation (4.4-20) shows that the wave shape is proportional 
to the value of the normal derivative of dU/dN of the trial mode U evaluated at the free surface, and 
Eq. (4.4-33), shows that the choice of the mass-distribution function M(S) has a strong influence 
on dU/dN. Hence, the selection of M(S) is a critical part of the numerical method. 

There are several conditions that M(S) should satisfy to make the trial functions U a good 
approximation to the true potential d>. Since the wall is impermeable, dU/dN is identically zero on 
the tank wall; thus, if y e were also zero, continuity would require that dU/dN on the free surface at 
the contact line should also be zero.. For numerical work, y c cannot be chosen to be exactly zero 
but it can be chosen to be small. Consequently, dU/dN on the free surface at the contact line should 
be "small," and the first requirement on M(S) is that it have a small numerical value at 5 = S c . By 
symmetry, both U and the wave amplitude are zero at the centerline of the tank 5 = 0. Hence, the 
wave amplitude must have a maximum between the centerline and the wall of the tank. Thus, the 
second requirement on M(S) is that it must cause the wave to have its maximum amplitude 
somewhere between 5 = 0 and 5 = 5 e . A mass distribution that satisfies both requirements is: 

M(5) = l+e-(5/5 e f (4.4-36) 

The small parameter e determines how closely dU/dN approaches zero at the contact point, and the 
exponent m helps fix the overall shape of the wave. 

The parameters e and m must be interrelated in order to satisfy the contact line condition, 
Eq. (4.4-22). To derive this condition it is assumed that U is proportional to 5 for the fundamental 
mode (which is known to be approximately true for high-g sloshing and, as shown by the following 
numerical examples, is also true for low-g sloshing). With this assumption, Eq. (4.4-22) can be 
manipulated to give: 

siny. 

S c [(dy/dS ) cos y c - (d%/dE)] + sin y c 

The term in curly brackets can be computed directly from the shape of the free surface. Equation 
(4.4-37) is not correct for the higher order trial modes, since they are not even approximately linearly 
proportional to 5. Fortunately, the contact line is always satisfied reasonably well when y c is small, 
so the discrepancy between e and m for the higher modes is not a serious limitation. 


W 


(4.4-37) 
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4.4.10 Numerical Examples 

Several numerical examples were computed for a spherical tank to demonstrate the use of 
the analytical and numerical methods described above. To begin the calculations, a FORTRAN 
computer code was written to predict the equilibrium free surface shape F(S) by numerically 
integrating Eqs. (4.4-12) and (4.4-13), as a function of Bo, contact angle, and Fill percentage. This 
computer code is described in Appendix B. Predictions were made for Bond numbers of 1 and 2, 
fill levels of 25%, 50%, and 78%, and a contact angle of 5*. Figure 4.4-3 shows the predicted 
equilibrium surfaces. ( y c > 0 was chosen because a contact angle of exactly 0* cannot be treated in 
the finite element structural code. The reason for not choosing Bo = 0 as an example case will be 
discussed later. The 78% fill level was chosen, rather than the obvious choice of 75%, because for 
Bo = 2 the lowest part of the surface coincides with the center of the tank, r = 0, z = 0.) 

A pre-processor (GIFTS) was used to compute the nodal coordinates for the finite elements 
of the structural simulation, using as input the tank shape and the free surface shape F(S) predicted 
by the FORTRAN code. ADINA, a commercially available structural code, was used for the 
example simulations, but other codes, such as NASTRAN, would have been equally suitable. The 
finite element model is shown in Figure 4.4-4. Because of the cos 0 dependency of the desired trial 
modes, only one-quarter of the tank had to be simulated, with the lateral surface aligned with 0 = 0° 
being a no-flow boundary dU/dN = 0, and the lateral surface aligned with 0 = 90° being an 
anti-symmetrical boundary (7 = 0. It was found from preliminary numerical examples that 20 - 25 
nodal points along the 0 = 0 line were more than sufficient to ensure convergence. 

All the wave shapes of the predicted trial modes were visualized with the aid of a graphical 
post-processor in order to select the desired subset of cos0 modes from the complete set of 
cos(2A + 1)0 modes which satisfy the imposed conditions for the one-quarter tank model. 
Generally, three such modes were selected. A QUICK BASIC computer code was written to perform 
the numerical integration of Eq. (4.4-35) and to compute the actual sloshing modes and the 
parameters of an equivalent mechanical model; this code is also described in Appendix B. 

Numerical experimentation was required to find suitable values for the parameter m. Wave 
shapes predicted by previous finite difference analyses [20] allowed a reasonable starting value to 
be chosen for some cases, such as m = 2.5 (and e = 0.1) for a 50% full tank with Bo = 1. An iterative 
process was used to refine the starting values; the first set of trial functions was used to determine 
an estimate of the true wave shape, from which a second estimate of m and e was obtained, and so 
on. In practice, only one or two iterations were required. Perhaps an even better procedure would 
have been to use the first estimate of the wave shape to determine a new, discrete distribution of 
M(S) for the second and later iterations, rather than continue to use Eq. (4.4-36); however, this 
procedure was not used for these example numerical calculations. 
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Figure 4.4-5 shows typical results for the trial mode predictions, in this case for Bo = 1 and 
50% filling level. The corresponding trial eigenvalues Q 2 were 0.4565, 1.5375, and 2.8853. The 
fundamental low-g slosh eigenvalue for this case was computed to be Q 2 = 0.7064 and the 
corresponding eigenvectors were b 2 lb, - 0.0158 and b 3 /b, = -0.0028. The eigenvalue of the second 
slosh mode was computed similarly to be = 14.306 and the eigenvectors were b,lb 2 = 0.2611 
and b 3 lb 2 = -0.0147. Figure 4.4-6 shows the computed slosh wave heights Hi and H 2 for this case, 
along with the wave heights (dUydNy&i and (dUydNyQ? of the first two trial functions, all 
normalized to have a peak amplitude of one. 

TABLE 4.4-1. SUMMARY OF LOW-G SLOSH PARAMETERS 
FOR A SPHERICAL TANK 



Fill % 

co 2 /(l +Bo)o/pR 2 

c» 2 /(g/B 0 ) 


U2R 0 

k s /oR 2 


25 

0.667 

1.335 

0.210 

0.749 

0.330 

Bo = 1 

50 

0.706 

1.413 

0.200 

0.708 

0.593 


78 

1.013 

2.026 

0.130 

0.494 

0.419 


25 

0.738 

1.106 

0.308 

0.678 

0.437 

Bo =2 

50 

0.816 

1.228 

0.250 

0.613 

0.642 


78 

1.221 

1.832 

0.168 

0.410 

0.450 


25 

1.299 

1.299 

0.745 

0.385 

0.601 

Bo = oo 

50 

1.573 

1.573 

0.580 

0.318 

0.772 


78 

2.193 

2.193 

0.350 

0.228 

0.521 


Note: = ^np R] x [filling fraction] 


Table 4.4-1 summarizes the computations of the fundamental slosh frequency and the 
mechanical model parameters, together with previous results [1] for high-g sloshing (i.e., for Bo = 
«*»). For comparison purposes, the predicted low-g frequencies are also presented in the conventional 
high-g non-dimensional form (£> 2 /(g/R 0 ). As can be seen, the non-dimensional eigenvalue 

Q 2 = co 2 /(l +Bo)a/pR \ for a given fill level decreases as Bo decreases. From physical reasoning, 
it is expected that the eigenvalue would be exactly zero for Bo = 0 when y c = 0; for these examples, 
the contact angle was 5*, so Q. 2 is probably not quite zero for Bo = 0, but the excessive amount of 
finite element computations required for Bo = 0, for which the static free surface is nearly a total 
spherical bubble, prevented this conclusion from being verified numerically. It should be noted 
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that c d 2 /(g/R g ) does not decrease monotonically as Bo decreases but first decreases and then shows 
a slight increase for Bo = 1, as a result of the changing interaction of the stiffening effect of surface 
tension and the softening effect of the increased wave length, compared to a high-g slosh wave. 

The predicted low-g slosh masses shown in the table are much smaller than the 
corresponding high-g masses and decrease as Bo decreases as a result of the increasing influence 
of the term proportional to H c in Eq. (4.4-28). The pendulum length increases significantly as Bo 
decreases. Since the pendulum attachment point is the center of the tank (for the physical reason 
that liquid cannot be set in motion in a spherical tank by changing only its angular orientation), the 
line of action of the pendulum mass can be outside the tank when Bo is small; that is, l, > R g . 

Experimental results for low-g sloshing in spherical tanks [6, 27] are limited to the range 
Bo > 10 and y c = 0, so a direct comparison cannot be made with these predictions. The trend of the 
data indicates that (£> 2 /{g/R 0 ) decreases as Bo decreases, which is thus in agreement with the trend 
shown in Table 4.4-1 for Bo > 1. The trend of the experimental slosh masses is also in agreement 
with the predictions. Furthermore, the frequencies given in Table 4.4- 1 agree fairly well with those 
predicted by a previous finite-difference numerical analysis [20], although the results in [20] are 
consistently slightly smaller, for example, for Bo = 1, the predicted Q 2 given in [20] is 0.619, 
compared to 0.794 by the present method. The slosh masses given in [20] are negative for Bo less 
than about five, because of an error in the slosh force analysis, and thus cannot be compared to the 
present results. 

4.5 Spherical Pendulum Rotary Slosh Model 

4.5.1 Background 

It has long been recognized that for liquids in symmetrical tanks, there is a strong tendency 
for rotational motion to occur throughout the steady acceleration range, even though the excitation 
to the tank may be planar. Although confirmation of this tendency under controlled experiments 
for low gravity in orbit remains yet to be accomplished, Peterson [28] has indicated its presence 
even in suborbital experiments conducted in parabolic flight trajectories. Thus, analytical modeling 
of the rotary slosh problem and its partial verification through earthbound experiments is an essential 
prelude to design of orbital experiments whose objective is the final confirmation of this complex 
fluid behavior. 

Development of a prediction for rotary slosh via hydrodynamic theory has not yet been 
accomplished. However, insight into the problem has been sought through studies of a spherical 
pendulum, which has been recognized as a potential analog for rotary liquid slosh. Nevertheless, 
even this approach quickly results in a relatively complex dynamics problem. Significant analytical 
study of the spherical pendulum has been reported by Miles [29]. Some experiments reported by 
Tritton [30] show that this classical system can display a variety of motions, including chaotic 
responses for certain conditions. On the other hand, Kana [31] has shown that the most significant 
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amplitude rotary liquid slosh in a scale model Centaur tank is at least periodic in nature, and can 
be modeled by a compound system which contains both a spherical and a linear pendulum. This 
report presents an extension of the latter work, and includes improvements to a harmonic balance 
model, which can be used to develop approximate dynamic parameters for both the spherical and 
the linear pendulums, and to a numerical approach which can be used to refine these results to more 
exact values. Finally, the effects of low gravity on this type of response are also explored. 

The physical configuration for the compound pendulum rotary slosh model originally 
reported by Kana [31], is repeated for convenience herein in Figure 4.5-1. The major developments 
of this report deal with the spherical pendulum, although the results will affect parameters for the 
linear pendulum part of the model as well. A summary of pertinent expressions which allow 
determination of parameters for the spherical pendulum will first be given. 

The general dynamic equations for the spherical pendulum were derived as: 


0+ co^ — <]> cos 0 sin0 + 2co„£ e 0 + — cos<}>cos0 = 0 

k / * 


and 


(4.5- la) 


<|) sin0 + 2<])0cos0 + 2co /1 ^4>sin0 -ysin<{> = 0 
wherein there has been included: 

, s IW, , C. r c, 
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(4.5- lb) 


(4.5- lc) 


Expressions for the cross-axis effective weight W e (co) and the in-line effective weight W,(co) 
were derived as: 


W^colcoscor = - 


P ,WV 
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.2 


4> sin © sin 4> - 2co0 cos 0 cos \|/ 


and 


- 0 cos 0 sin y 




(4.5-2a) 


p,uy 

^(co) cos cor = i - LJ r 
x 0 (a 2 

-i-p.^coscor (4.5-2b) 

A harmonic balance approximate solution to Eqs. (4.5-la, 4.5-lb, and 4.5- lc) was first 
developed. Radial force equilibrium resulted in the expression: 


<J> sin0cos<{> + 2co0 cos Osin >|/- 0cos0cos\|r 
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tan 0 O - a sin 0 O = — 
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Circumferential force equilibrium resulted in the expression: 

C* = ±1 


^octo^ sin 4> 0 CCW 


V4 1 J 


sin0 o CW 


(4.5-3a) 


(4.5-3b) 


(4.5-4) 


A corresponding harmonic balance expression for the complex cross-axis weight was 
developed from Eqs. (4.5-2a and 4.5-2b) to result in: 


.2 


W,(co) = -p 1 W 1 /^ r sin0e‘* + 2cop 1 iy i /^cos0e'f y+ ^ 


+ P,W 1 /®cos0e ,> +p 1 U' ) * 


(4.5-5a) 


in which there was included: 

\jr = cor - (y 0 + 4> 0 ) = <J) — Yo 


(4.5-5b) 


and 


tanyo = —5, Cot = ~ 


l-a z ' "" 2 

From this, there was defined along the real cross-axis a co-phase component W cc as: 


(4.5-5c) 


W cc = P,W, 


/ , 

— sin 0 O sin 4> 0 + cos Q 0 F l 
x 0 


\(a) sin Xq = W £ ( co) 


and along the in-line axis a quadrature-phase component W CQ as: 


W CCI = -P,IV, 


— sin 0 O cos <J> 0 + cos 2 ©oF/a) cos 

*0 


* =P,W 1 -^(co) 


where W CQ is taken as positive to the left. In these expressions, there was also used: 

Xg = 2(y 0 + 4> 0 ) + 60 

and 


F,( a) = 


a 


(4.5-6a) 


(4.5-6b) 
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(4.5-7b) 
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The distortion angle Eq was included to allow the pendulum weight ratio p, to remain in 
the realistic range between zero and one. 

In [31], development of model parameters from the above equations was accomplished by 
the following steps: 

1) Circumferential damping with corresponding position lag angle <j) 0 and steady 

deflection 0 O were obtained by simultaneous iteration on Eqs. (4.5-3 and 4.5-4). The 
condition was that <t> 0 match experimentally observed values, where they could be 
discerned (i.e., fora > 1.0). Fora< 1.0, only guesses on <J> 0 could be made, since 0 O was 
too small to isolate. 

2) These values of t*, <]> 0 , and 0 O were then used in Eqs. (4.5-5, 4.5-6, and 4.5-7) to calculate 

cross-axis weight components. For this, an assumed value of radial damping and 
error angle Eq were varied until the calculated weight components matched those which 
were measured for a given excitation frequency during the slosh experiments, and the 
resulting weight ratio P, was a realistic value. 

With the above steps, all parameters for the spherical pendulum were first determined from 
the approximate solution. Thereafter, the damping values and natural frequency G)„ were used 
in Eqs. (4.5- la, 4.5- lb, and 4.5-lc) for a numerical time-step solution for 0(r) and <J>(r). These 
solutions were then subsequently input to Eqs. (4.5-2a and 4.5-2b) to obtain more accurate weight 
components as functions of time. However, in doing this, the results were compared on a magnitude 
basis only, with phase being ignored. 

The primary objective of the present work is to develop more accurate parameters for the 
spherical pendulum part of the slosh model by use of similar numerical solutions of the governing 
equations. However, both magnitude and phase of the cross-axis weight components will be 
included. As a result, parameters for the linear pendulum part of the compound slosh model will 
also be shown to be affected. Thus, a more accurate slosh model which matches all dynamic 
properties of the cross-axis and in-line weights will result. 

4.5.2 Revised Approximate Model 

As indicated in the previous section, judicious choices of damping parameters must first be 
made for the spherical pendulum before a direct numerical solution of the governing equations can 
be attempted. These values must be such that the solutions for 0(f) and 4>(f) will produce cross-axis 
weights which can be made to match the experimental slosh data when a plausible mass ratio (3, is 
included. Therefore, to assure that judicious initial values are selected, the harmonic balance 
approximate solution is still first used to estimate values of damping, deflection, and position lag 
angles. However, in the present work, the model is revised so that £o is set to zero in Eq. (4.5-7a). 
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By further search, it was found that plausible solutions for p, indeed could be obtained by significant 
revision of the damping values i * and £e,. Otherwise, the approximate parameters are all determined 
using the same detailed approach described previously in Steps 1 and 2. 

The results obtained from the revised approximate model with £<, = 0 are shown in Figures 

4.5-2 - 4.5-5. Note that numerical solution results are also given, but will be ignored until later. 
Here the approximate results are presented as continuous solid and dashed lines, although it is 
understood that values were obtained only at frequencies where experimental data were available. 
By comparing these results with those presented in [31], it can be seen that for values of a < 1.0, 
significantly greater values of damping for both and ^ are necessary to develop a plausible 
model. At the same time, the steady deflection angle 0 O was not changed very much, but the values 
of position lag angle <{) 0 and mass fraction p, change significantly. Even so, they remain in a plausible 
range (i.e., P, < 1.0). The net result of this development is that a better model of the experimental 
slosh data can be achieved, as will be described further with the numerical approach. 

4.5.3 Numerical Model 

As previously indicated, the objective of the numerical model is to seek solutions for 0(f) and 
<|>(f) so that matching of all dynamic properties of the experimental slosh data can be achieved. 
Since the spherical pendulum alone produces the cross-axis weights, it is given attention first, and 
then the linear pendulum is established by means of the in-line weights. However, in developing 
the approach to the spherical pendulum, it is first appropriate to provide some preliminary discussion 
about the expected forms of the solutions and to modify the cross-axis weight expressions. 

Forms of Solutions 

Steady state polar solutions for 0(f) and <t>(f) from Eqs. (4.5-la and 4.5-lb) are of elliptical 
shape similar to the example shown in Figure 4.5-6a. This polar plot represents counter-clockwise 
motion, as will all the solutions considered in this report. (Discussion of clockwise solutions will 
be given in the conclusions). Furthermore, a modified notation (compared to that used in [31]) has 
been employed in Figures 4.5-6a and 4.5-6b in order to allow more generalization, and to emphasize 
in which plane a given variable is defined. 

Solutions of interest are periodic in cof , and in order to study the details of a single period one 
can pick an initial reference time: 

f 0 = 2nn/(a (4.5-8) 

where n is an integer. By referring to Figures 4.5-6a and 4.5-6b, we can define: 

4>,(f) = cof-<}> 0 (f) (4.5-9a) 

and 

= “f — <t>oi(f) (4.5-9b) 
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where 


(4.5-9c) 


4>oi(0 — OoCO 4*qj 

Furthermore, combining Eqs. (4.5-9b) and (4.5-9c) there results: 

4>x(0 = co[r - <j 0 (r)/co] - (4.5-10) 

From these equations, we further note that at t = t 0 \ 

4>o( r o) = 0 and 0 O i (r 0 ) = (4.5-11) 

Note that if the solutions of Eqs. (4.5-la and 4.5-lb) are such that the polar plot in Figure 
4.5-6a becomes circular, then Eq. (4.5-10) takes on the more familiar monotonically-increasing 
straight line form, in which 4 > 0 (r) is zero. However, if the solutions form an elliptical polar plot as 
shown in Figure 4.5-6a, it will be shown later that <)> 0 (r) is a periodic function superposed on the 
otherwise straight line <j \(t) = (0 1 — (j)^. 

Modification of Cross-Axis Weight Expressions 

It is now appropriate to use Eqs. (4.5-2a and 4.5-2b) to develop effective weights as 
functions of time so that comparisons with experimental data can be made. However, as given, 
these expressions include some restrictions which were useful for the approximate model solution, 
and will be eliminated for the numerical model. In particular, the last two terms in the brackets on 
the right-hand side of Eqs. (4.5-2a and 4.5-2b) were defined in terms of a rotating pendulum whose 
plane of oscillation lags at a constant angle relative to the space-plane angular displacement 
<(),(r). As indicated in Figures (4.5-6a and 4.5-6b), it was assumed that: 

Ya. = MO = Yo-<l>a, (4.5-12) 


where <j> x (r m ) is the angle between the apex of the ellipse and the x-axis. To account for coupled 
radial [i.e., 0(r)] and rotational [<K01 effects, the angle y ( (0 was defined by Eq. (4.5-5b) as: 

V,(0 = cor - hfo + 4>o(01 (4.5-13) 


and using Eqs. (4.5-9a, 4.5-9b, and 4.5-9c) there results: 

V,(0 = WO-Yo + ^t, (4.5-14) 

Use of this approach only partially accounts for coupling effects that are present, and will now be 
generalized for development of the numerical model. 

Equations (4.5-2a and 4.5-2b) are considered directly for a pendulum which is rotating 
with variable angular velocity <j>(r) and with a plane of oscillation which makes the angle 4>,(r) 
relative to the excitation plane. For this, Eqs. (4.5-2a and 4.5-2b) become: 

— ftwy 

1^(0 cos cor = — 

x 0 a> 


4> sin 0 sin (])- 2<})0 cos 0 cost])— 0cos0sin0 


(4.5-15a) 


and 
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W,(t) cos tor = 


jW 

XqOl ) 2 


.2 


4> sin0cos(f) + 24)0cos0sin<J)— 0cos0cos<}) 


+ p 1 W 1 coscof (4.5- 15b) 

In which it is understood that <j) = §,(t) and 0 = 0,(r) are each calculated from coupled numerical 

time-step solutions of Eqs. (4.5-la and 4.5-lb), and W e (t) and W,(t) are vectors which have a 
time-phase relative to cos cor. 

Thus, we have: 

iv e (r)coscor = 1 W^OIe'^coscor (4.5- 16a) 

and 

W,(r)coscor = |M 7 t (r)|e" 75 °'coscor (4.5- 16b) 

where ^ and ^ are lag angles for the respective weight components. These lag angles are found 
by: 

e = ©Oc-'o) (4.5-17a) 

$ 0 . = co(r < -r 0 ) = ^ 0e + 7t/2 (4.5-17b) 

where t 0 is given by Eq. (4.5-8) and t c and t, are the times for a maximum occurrence of the respective 
weight just succeeding the reference time t 0 . (Note that if a preceding maximum occurrence is used, 
then ^ becomes a lead angle rather than a lag angle). 

When the cross-axis weight solutions are developed according to Eqs. (4.5- 15a and 
4.5- 15b), it is found that a time lag error exists relative to the experimental weight measurements. 
This error is found to be dependent on the variation of <j) 0 (r) with time. In view of Eq. (4.5-10) we 
note that: 

W £ [<t>,(r)] = W e {<j>,[*-<j> 0 (r)/co]} 

Guided by this, we set: 

t' 0 = t 0 -<p 0 (t c )/a (4.5- 18a) 

This results in a lag angle: 

5 o c = ©& " O = $oc + UQ (4.5- 18b) 

Therefore, since the cross-axis weight is expressed as: 
w c (t) = W CC -JW CQ 
there results: 

W cc = \W c (t)\ cos^ and W CQ = \W e (t)\ 005(^-71/2) (4.5-19) 
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Similarly for the in-line weight: 

W t (t) = W. R -jW t , 
there results: 

W.* = \WM cos ^ and W tl = \W t (t)\ cos(^-7t/2) (4.5-20a) 

where 

S'o. = - 1 \ 0 ) = ^ + MO = t'oc + (4.5-20b) 

Spherical Pendulum Procedure 

The procedure for development of model parameters will be described in steps along with 
results for a specific frequency point. For illustration, the point at a = 0.9720 was selected, which 
includes the following experimental slosh data (see [31]): 

W cc = 108.0 lb(480.4 N) and W CQ = 59.57 lb(265.0N) 

1) For this case. Step 1 of the approximate model provides values of damping as 
C* = 0.017 and Cei = 0-0 22 

These values appear as part of the solid and dashed lines, respectively, in Figure 

4.5- 2. However, after trial in the rest of the steps to follow, it was found that even 
more damping was necessary to produce sufficient cross-axis force. Therefore, these 
values were increased to: 

C* = 0.025 and Ce, = 0.045 

so that ^ = 0.090 

2) The above values of damping were used along with the DYSIM (Dynamic 
Simulation) computer program to compute 0(r) and <J>(f) from Eqs. (4.5- la and 

4.5- 1 b). This program is based on a fourth-order Runge-Kutta numerical integration 
scheme. The results are plotted for reference in Figures 4.5-7, 4.5-8, and 4.5-9 for 
the time period 45 to 50 seconds. This time was sufficient to establish steady-state 
after using the following initial conditions: 

0(0) = 0.0038 radian <J>(0) = —1.534 radian 

0(0) = 0.000 rad/sec (J>(0) = 2.7300 rad/sec 

In this and all cases developed herein, the time step increment was 0.002 sec. 
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For convenience 4> 01 (O was also computed from Eq. (4.5-9b) as: 

4> 01 (0 = cor - ^(r) 

and plotted in Figure 4.5-1 0. It is clearly seen to be a periodic function for this case. 

4) From the above results, the following parameters are obtained from Eqs. (4.5-8 and 
4.5-11) for n = 21 periods: 

r 0 = 48.331 sec ^ = 0.42 rad 

5) The above results for 0(r) and <K0 are now used to compute W c (t) cos (tit and 

W, (r) cos cor from Eqs. (4.5- 15a and 4.5- 15b). The results for the periodic component 
filtered at frequency to are plotted in Figures 4.5-11 and 4.5-12, respectively. For 
this computation, the value of p t is adjusted so that the magnitude of | W c (r)| equals 
that for the experimental data: 

l 

\W e (t)\ = [W 2 cc + W 2 Q f = 123.3 lb (548.4 N) 

This results in p, = 0.160. Furthermore, from the results in Figure 4.5-1 1 we obtain 

t c = 48.162 sec. With this and the results from Figure 4.5-10 used in Eq. (4.5-9c), 
we obtain <J) 0 (r c ) = 0.94 radian. Therefore, from Eq. (4.5- 18a) we obtain r' 0 = 47.987 
sec, and from Eq. (4.5- 18b) there results ^ / 0c = 27.4*, which agrees with the 
experimental value. 

Thus, a very close matching of both magnitude and phase has been achieved for the 
cross-axis weight data, and has resulted in a very plausible value for the weight ratio 
P,. Generally, this is accomplished only after some additional adjustment of the 
damping ratios, as was indicated above in Step 1. 

6) The in-line components for the spherical pendulum are obtained from W t (t) cos (tit, 

which is plotted in Figure 4.5-12. Note that the scale on this plot has been adjusted 
according to P, = 0.160, as found in Step 5 above, so that Figure 4.5-12 also provides 
|W,(f)| =107.7 lb (479.1 N). The timer, =48.70 also is established from this plot. 
Therefore, from Eq. (4.5-20b) the in-line lag angle is obtained as ^ = 1 1 1.8*. From 
this and Eq. (4.5-20a), there results W,*=-40 lb (177.9 N) and W tl = 99.5 lb 
(442.6 N). 

This information is required for determination of the linear pendulum parameters, 
as will be shown later. 
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Spherical Pendulum Results 

Spherical pendulum model parameters were similarly calculated for several other frequency 
points, and the results plotted as discrete points in Figures 4.5-2 - 4.5-5. Generally, it can be seen 
that the resulting damping ratios ^ and £e, and mass ratio p, all are higher than the corresponding 
values estimated from the approximate model, only for the frequency range a < 1.0. This could be 
expected, since for a >1.0 the polar plots rapidly become circular, and the approximate and 
numerical solutions become more identical. Even so, the approximate values are useful in all cases 
to use as initial values for numerical model development. Note also that the steady deflection angle 
0 O is essentially the same for both cases (for the numerical model 0 O is taken as the average value 
from the 0(f) - plot). Finally, from Figure 4.5-4, it can be seen that the values for position angle <|> 0 
also deviate somewhat only for a < 1.0. In this case, <t> 0 is taken as: 

4*0 = $0j 

i.e., the value for <[> OI (0 at t = t 0 . For a > 1.0 and more circular orbits, <J) 0 i(0 becomes constant and 
<P(r) in Figure 4.5-9 becomes a straight line. 

Some further examples of types of polar plots are shown in Figures 4.5- 1 3 - 4.5-15. Figures 
4.5-13 and 4.5-14 show samples of plots from actually developed model data given in Figures 4.5-2 
- 4.5-5. From Figures 4.5-7, 4.5-13, and 4.5-14, it can be seen that dramatic changes in the orbits 
occur (i.e., in magnitudes and shapes), which correspond qualitatively with what is observed 
experimentally in the liquid behavior. These changes result from variation of damping values as 
well as frequency. For example, Figure 4.5-15 shows several orbits computed for fixed excitation 
frequency a = 0.972 and circumferential damping ^ = 0.017, but with varying radial damping £ei. 
For large £e,, the orbit is nearly circular, so that a relatively large cross-axis weight would be 
produced. For small £e,, the orbit reduces to the limiting case of a linear pendulum, and zero 
cross-axis weight results. For this also <}>(t 0 ) becomes essentially zero, so that the in-line weight 
reduces to that for a linear pendulum with amplitude and phase determined by £e,. Furthermore, 
the function <J >(r) approaches a straight line for the more circular orbit, while corresponding results 
for the linear pendulum form a stair-step function. In the other cases, the results are similar to those 
of Figures 4.5-7 - 4.5-10. 

Linear Pendulum Discussion 

In [31], the parameters for the linear pendulum in the compound model are developed 
totally from the experimental data. Herein, this approach is changed to recognize that for less 
circular (i.e., more elliptical) orbits, the corresponding cross-axis and in-line components of the 
spherical pendulum are not equal, since (J) is not a constant for such orbits. That is W cc * W t , and 
W CQ * W tR . This conclusion is evident from the results of the numerical model. Therefore, the 
combined system weight equations are now formed as follows: 
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(4.5-2 la) 
(4.5-21b) 


W K = W tR + + P 2 W^jF L (a)cosy 0 

W, = ^ / + p 2 M',F t (a)sin Vo 
where the total complex in-line weight is: 

W, = Wg-iW, 


a 

F l { a) = , 

(4.5-22a) 

[( l-a 2 ) 2 + 4&a 2 f 


2 

(4.5-22b) 

tan Yo2 = : 2 

1 - a 


Equations (4.5-2 la and 4.5-2 lb) indicate that the total in-line weight is produced by the 
two pendulums. The contribution from the spherical pendulum results from the numerical solutions 
for W eR and W d as indicated previously for the spherical pendulum procedure. The contribution of 
the linear pendulum results from its amplitude given by Eq. (4.5-22a) and phase given by the angle 
\|/ 0 . However, if there is some coupling between the two pendulums required for matching the 
in-line slosh data, then there will result: 

Y 02 * Vo 

For this case, \|/ 0 is obtained from combining Eqs. (4.5-2 la and 4.5-22b) to obtain: 


tan\)/ 0 


w r -w, r -&w x 


Then a coupling difference angle Eo 2 will result such that: 


(4.5-23) 


£o 2 = V 0 -Y 02 (4.5-24) 

Linear Pendulum Procedure 

Parameters for the linear pendulum can now be developed by means of the above results. 
For illustration, the data for the frequency point a = 0.972 will be continued. For this, the total 
measured in-line weights (see [31]) are: 

W R = 389.61b (173.9 W) and W, = 343.31b (1527.0 N) 

1) Given the above total measured in-line weight values W R and W, and the in-line 
weight components W tR and W eI computed from the numerical model one may then 
assume a value for For this case, this damping value will be taken as constant 
for all frequencies at: 
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Ce2 = 0.010 


With this, F l ( a) is computed from Eq. (4.5-22a) and y 02 is computed from Eq. 

(4.5-22b). For this example: 

F l ( a) = 16.14 and y 02 = 19.40° 


2) Equations (4.5-2 la, 4.5-2 lb) are now manipulated to where is eliminated from 


them. This results in: 
(W r -W. r ) 


(W,-W tl ) 


sin\j/ 0 -cos\)/ 0 


1 


F l { a) 


(4.5-25) 


With the above data, this equation is solved for \|/ 0 by trial and error. For the example 

case, there results: 

% = 30.7° 


and from Eq. (4.5-24) this results in: 

£o 2 = 11.3° 

Thus, a non-zero angle indicates that some degree of coupling between the two 
pendulums is necessary to match the in-line weight data. Furthermore, the linear 
pendulum mass ratio can then be obtained from Eqs. (4.5-21a) or (4.5-2 lb) as: 

P 2 = 0.310 

Linear Pendulum Results 

Further results for the linear pendulum which correspond to the data developed for the 
numerical spherical pendulum model are shown in Figure 4.5-16. The resulting mass ratio 0 2 varies 
more or less similarly to p, for the spherical pendulum. The coupling angle £o 2 varies significantly 
over the frequency range, which indicates a corresponding large variation of coupling between the 
two pendulums, as long as is held constant. 

It is appropriate to raise the question of whether this coupling can be reduced to zero 
providing that is allowed to vary. It was found that this approach led to impossible values of 
damping at the lower frequencies. Therefore, for simplicity, the ^ constant value was selected. 
This may or may not be the optimum approach for every set of experimental data developed. 

4.5.4 Effects of Low Gravity 

The preceding developments have concentrated on an approach for developing parameters 
for a compound slosh model that can match the effective weights measured for liquid motions in a 
tank which is subject to a given steady acceleration. The experimental data was acquired in an 
earthbound system, and the equations were derived for a spherical pendulum in an earthbound 
system. The next logical question deals with whether the model developed can be used to predict 
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results that might be expected in a low gravity environment in orbit The methods developed for 
both the harmonic balance approximate model and the numerical model will be used to shed some 
light on this question, by showing how low gravity can affect the response of the compound slosh 
model. 

An immediate effect of low gravity on both a spherical and a linear pendulum is that the 
natural frequency is dramatically reduced. For a liquid, this effect is so pronounced that surface 
tension forces then dominate, rather than gravity forces, as has been described earlier in this report. 
For a pendulum in low gravity, this corresponds to adding a weak, torsional spring at its support. 
Thus, if no other forms of nonlinearity enter the problem and if one can assume that excitation will 
still occur in the vicinity of the natural frequency, however low it may be (i.e., near a = 1.0), then 
the essential modeling approach developed herein remains applicable. However, a very important 
effect on the results will additionally occur because of changes in damping, as the previous 
development shows that dramatic differences in response of a spherical pendulum occur for a given 
frequency when either ^ or ^ change by only slight amounts. Although the effects of low gravity 
on damping in a given tank are not yet completely understood, it is known that damping tends to 
increase as gravity decreases [i.e., Eq. (1.1-1)]. Therefore, the subsequent discussion concentrates 
on what can happen to a spherical pendulum excited near a = 1.0, when the damping and are 
allowed to increase. 

Equations (4.5-3, 4.5-4) which represent part of the harmonic balance approximate model 
are first used for this purpose. Equation (4.5-3a) was first solved for a variety of position angles 
<{> 0 for incremented values of frequency ratio near a = 1.0, and the corresponding values of steady 
deflection 0 O were noted. This results in the plot shown in Figure 4.5-17. Then, corresponding 
values of were calculated from Eq. (4.5-4), and the results plotted in Figure 4.5-18. By using 
the two figures together, one can estimate what type of responses occur at a given frequency for an 
initial value of damping and also how the character of the response changes as this value of 
damping is increased. 

As the previous model developments indicate, generally an increase of damping ^ leads 

to larger values of <]) 0 , and therefore, correspondingly larger values of cross-axis effective weights, 
all of which is independent of £e, only for the approximate model. However, this alone is not the 
whole story. As was previously shown, at a fixed value of C* the character of the response also 
changes with £ 0 , increase. Therefore, use of the numerical model is ultimately necessary for a more 
exact determination of response changes that occur due to variations in damping that result from 
low gravity. 
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Relating these results to changes in gravity may be accomplished by estimating the damping 
change according to Eq. (1.1-1). However, the computed value is most directly associated with the 
radial damping £e, in this model. At present, it must be assumed that some similar relationship 
exists for the circumferential damping C*. On the other hand, the problem is even more complicated 
than such an approach may imply. By programming various combinations of damping into Eqs. 
(4.5-la and 4.5-lb), for the spherical pendulum, it was quickly found from numerical solutions that 
steady state response occurs only within certain regions of damping — otherwise drifting or chaotic 
responses occur, as predicted by Tritton [30]. For the present case, it was found that stable periodic 
motions could be found for the ranges of frequency, damping and associated phase angles indicated 
in Figures 4.5-17 and 4.5-18. However, these values were found only by trial and error, and more 
work is required before a better understanding of the behavior is possible. Furthermore, 
establishment of more accurate relationships for both and £ei as functions of low-g are required 
before any predicted results can be meaningful. 
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5.0 GROUND-BASED SUPPORTING EXPERIMENTATION 

5.1 Introduction 

Several different studies were conducted to help establish requirements for flight 
instrumentation. The response characteristics of liquid-vapor interface sensors were investigated 
by laboratory experiments to determine an appropriate set of spatial and frequency response 
requirements; these experiments used non-flight hardware and instrumentation. Load cell designs 
were surveyed to determine if designs exist that can measure the small slosh forces anticipated for 
a flight experiment. This section summarizes those investigations. 

5.2 Ground Tests Of Liquid-Vapor Sensors 

5.2.1 Background 

The geometric description of the static and dynamic free surface is among the most useful 
fundamental information that can be obtained about liquid motions in moving tanks. This is 
especially true for low-gravity conditions, since the shape of the free surface slosh wave is the only 
evidence obtainable from tests that is directly related to the contact line condition (i.e., to the dynamic 
contact angle). In addition, the natural frequency of the motions can also be determined from the 
time history of the dynamic free surface location, and the damping can be obtained by the time-decay 
of the free surface wave motion after the tank motion itself ceases. All this information obtained 
from surface location measurements is needed to validate and improve analytical models and to 
acquire fundamental understanding about the relevant surface physics. Thus, it is imperative that 
the static and dynamic free surfaces be "visualized" in some way during flight experiments. 

The most common visualization method used in laboratory studies of sloshing is a 
combination of; (1) cinema or video recording of the free surface motion; and (2) probes that 
measure the time history of the free surface location at one or more points. Cinema or video recording 
is an appropriate method of obtaining a qualitative overview of the motion and of selecting conditions 
for further study, while an array of surface probes is an appropriate method of obtaining quantitative 
data about frequencies, damping, resonant conditions, etc. For the COLD-SAT flight experiment 
defined in Section 3.2, cinema or video recording from outside the tanks is not possible because 
the COLD-SAT cryogenic tanks are not transparent; internally-mounted cameras are also not 
possible because of the substantial insulation required for the camera and the need for a source of 
lighting. For the Shuttle-based experiment defined in Section 3.3, the use of external cameras is 
practical and is therefore included in the plan. For both experiments, liquid-vapor interface probes 
are proposed as a method to obtain quantitative data. 

5.2.2 Liquid-Vapor Sensors 

Liquid surface sensors of the type commonly used in laboratory studies of sloshing cannot 
be used in the "weightlessness" of low gravity because they depend on the weight of any liquid 
adhering to the sensing elements to remove the liquid rapidly when the probes exit the liquid. 
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Instead, some other mechanism must be used to remove the adhering liquid so that a time-correlated 
signal can be generated which indicates when the sensing element is no longer immersed in liquid. 
At least one space experiment has used a temperature sensing element containing an electrical 
resistor for this purpose [32]. When the sensor was immersed in liquid, the heat generated by the 
electrical current in the resistor was conducted away more rapidly than when it was immersed in 
vapor, consequently, the temperature of the sensor was lower in liquid than in vapor, and the change 
of temperature thus indirectly indicated when the element entered and exited the liquid. These 
kinds of devices have been investigated more thoroughly for the COLD-SAT experiments [33, 34]. 
In addition, fiber optic devices have also been proposed as liquid-vapor surface sensors [34, 35]; 
Figure 5.2-1 shows the principle upon which these sensors operate. 

Although the temperature response time of most resistive sensors is slow (on the order of 
five to ten seconds), the required response time for the slosh flight experiments is not particularly 
demanding; Tables 3.2-2 and 3.3-2 indicated that the slosh waves for the two defined space 
experiments all have periods greater than 70 seconds. Consequently, the response time of resistive 
sensors should be adequate. In laboratory tests, the response of fiber optic sensors is practically 
instantaneous; it remains to be shown, however, whether liquid adhering to the probe tip in low 
gravity can degrade the performance of the probe. 

It is concluded that resistive sensors, and possibly fiber optic sensors, can be adapted to 
the requirements of determining the location of the liquid-vapor interface in low gravity flight 
experiments. The number and location of such probes required to obtain an accurate resolution of 
the shape and motion of the surface were therefore the objectives of the present laboratory 
experiments. 

5.2.3 Ground Tests of Liquid-Vapor Sensor Arrays 

An existing 1/5-scale model Centaur G-Prime tank was available for the laboratory studies. 
Water was the test liquid. Near-resonant planar sloshing was excited by oscillating the tank with 
a horizontal shaker. The test apparatus is sketched in Figure 5.2-2. 

Modified "wheatstone" wave height transducers [36] were used in the tests to simulate the 
resistive or fiber optic sensors of a flight experiment. Each such liquid-vapor sensor was made of 
a pair of thin insulated conductors, slightly separated, with a short section of each conductor exposed 
for contact to water, and connected electrically across part of the resistance in one leg of a wheatstone 
bridge. The short exposed section of the sensor was oriented horizontally (i.e., parallel to the static 
liquid surface) so that a distinct indication could be obtained of the time when the liquid surface 
passed the sensor. An individual probe was composed of four sensors mounted on a thin vertical 
rod and separated vertically by one inch (2.54 cm). Three probes (Rl, R2, and R3) were constructed 
and mounted in the tank on a brace along the tank diameter as shown in Figure 5.2-2. The entire 
array of probes could be adjusted angularly relative to the tank excitation direction, and vertically 
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to position the sensors at any desired location relative to the static liquid level. About seventy-five 
slosh tests (not counting repeats) were conducted for a number of different orientations of the probes 
and for three different slosh wave amplitudes. No attempt was made to simulate the curvature of 
the static liquid surfaces that will occur in low gravity. 

Figure 5.2-3 shows a reproduction of the strip-chart recording of the output of the probe 
array from a typical test. When the slosh wave passed upward across a sensor, the sensor went 
from a "dry" state to a "wet" state, and there was a consequent abrupt increase of the bridge output 
voltage. When the wave passed downward, the sensor went from a "wet” to a "dry" state, and there 
was a consequent decrease of the bridge voltage, which was, however, not as abrupt as for the 
upward passage since some water adhered to the sensor and later dripped off. The magnitude of 
these voltage jumps is related to the sensor sensitivity and calibration, but not to the height or depth 
of the wave above or below the sensor. 

Slosh frequency — The time period between successive dry or successive wet indications 
for a given sensor, which is, of course, the slosh wave period, was repeatable to high accuracy. As 
an example, the average period (derived from the known chart travel speed) for the test results 
shown in Figure 5.2-3 was 1.205 seconds, and the maximum variation from one cycle to another, 
or between sensors, was 0.005 seconds. The actual slosh period, which was set by the shake table 
frequency, was 1.205 seconds. It was concluded that digital sensors are easily capable of establishing 
slosh frequencies. 

Static liquid level and slosh wave shape — The output of the entire array was available 
to establish the static liquid level and the shape of the slosh wave. As an example of one method 
that can be used to interpret the sensor data to obtain this kind of information, the test that yielded 
the data shown in Figure 5.2-3 will be analyzed. Figure 5.2-4a shows the configuration of the probe 
array, which was aligned with the tank excitation direction. The static liquid level was halfway 
between the middle two sensors of each probe. The initial uncertainty in the static liquid level 
(assuming that there was no visual evidence, such as would be the case for COLD-SAT experiments) 
was, therefore, ±0.5 in (1.27 cm) because the level could have been anywhere between adjacent 
wet and dry sensors. 

Since the sensors at a given vertical level did not all indicate "dry to wet" or "wet to dry" 
at the same time, a technique was developed to interpolate the discrete data from each of the four 
sensors of a probe to determine a continuous time history. The slosh motion is periodic, so an 
appropriate interpolation technique was to fit a least-squared-error sine wave to the data: 

ri=A 0 +A 1 sin(2ro/T / + p 0 ) (5.2-1) 
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Here, rj is the position of the liquid surface relative to lowest sensor of the selected probe, x, is the 
period of the slosh motion (1 .205 sec for this example), and 14 is a phase shift relative to probe Rl. 
The parameter A„ is interpreted as the position of the static liquid level above the lowest sensor of 
the probe and A, as the amplitude of the slosh wave at the radial location of the probe. The data 
used to derive the least-squared-error sine wave consisted of: (1) the time (after an arbitrary starting 
point) when each sensor went from dry to wet and from wet to dry, and (2) the vertical position of 
each sensor relative to the lowest sensor. The time data were obtained from the strip chart recordings 
of each test by measuring the length intervals between voltage "jumps" which were then converted 
to seconds. As can be inferred from Figure 5.2-3, there was almost no ambiguity in the dry-to-wet 
measurements and very little in the wet-to-dry measurements; generally, measurements for three 
or four slosh periods were averaged to remove what ambiguity existed. The position data for the 
sensors were determined from direct measurements; as indicated earlier, the corresponding sensors 
of each probe were at the same vertical position, and each sensor was separated vertically by 1 in 
(2.54 cm). 

For this example, the results of the interpolation for the three probes are summarized as: 

Probe Rl: A 0 = 1.701 in A, = 1.859 in n 0 = (T 

Probe R2: A 0 = 1.511 in A, = 1.972 in *i o = -0.06’ 

Probe R3: A a = 1.440 in Aj = 1.974 in ^ o = -0.75* 

The interpolation method established the location of the static free surface (by averaging the A 0 
data) as 1.551 in (3.940 cm) above the lowest line of sensors, and reduced the uncertainty in that 
position from ±0.5 in (1.27 cm) to ±0.1 10 in (0.279 cm). The true position is 1.5 in (2.54 cm) to 
within the accuracy of the test setting. 

The sine waves fitted through the data for each probe are shown in Figures 5.2-5a - 5.2-5c. 
The predicted shape of the slosh wave, obtained from the amplitudes A 0 and A; of each sine wave, 
is shown in Figure 5.2-5d. Considering the small phase differences |i 0 between the probes, the 
phase has been neglected in this composite wave shape. The predicted slosh wave shape is a 
reasonable approximation of the actual wave shape, and in fact, if the static surface level had been 
predicted to be flat, the prediction would have been even more realistic. It is apparent that the 
sloshing was slightly nonlinear, which was also visually observed during the test The predicted 
value of * 2 in (5 cm) for the slosh wave amplitude at the wall was somewhat smaller than the true 
upward value of 2.65 in (6.7 cm) and even slightly smaller than the true downward value of = 2.2 in 
(5.588 cm). A Fourier series interpolation scheme could account for the nonlinearity, which would 
further improve the estimate of the static liquid position and the wave shape. 
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Figure 5.2-6 shows an example of the improved results that can be obtained by using more 
sensors per probe. The predictions shown in this figure were actually obtained by combining the 
data from two identical slosh tests, but in which for one the array was positioned with the lowest 
sensor 0.5 in (1.27 cm) below the static surface (i.e., three sensors were exposed) and for the other 
the array was positioned with the highest sensor 0.5 in (1.24 cm) below the static surface (i.e., three 
sensors were submerged). This data set thus simulates an array in which there are six sensors along 
each probe, equally spaced at 1 in (2.54 cm) intervals. (Two of the sensors for each test overlapped, 
which is the reason that the simulation does not represent eight sensors per probe.) The slosh 
frequency and wave amplitude were the same as for the test discussed previously. The static liquid 
level was now predicted to be 2.550 in (6.480 cm) above the line of lowest sensors, compared to 
2.5 in (6.35 cm) of the test, with an uncertainty of ±0.068 in (0.172 cm), and the slosh amplitude 
at the wall was predicted to be » 2.4 in (= 6.1 cm), compared to 2.65 in (6.7 cm) for the test. 

Other tests were conducted with the sensor array not aligned with the tank excitation 
direction; typically, angles of 15*, 30', 45*, and 60* with respect to the excitation were used. The 
direction of the excitation can be computed from the data from two or more such tests (or from two 
or more arrays). It is assumed that the maximum wave amplitude at a given radial position varies 
with the cosine of the angle between the excitation direction and line of the array. For example, 
using the data from two separate tests in which the array was aligned at 15* and at 45' to the excitation 
direction, the excitation direction, and thus the line of the peak wave amplitude, was predicted to 
be 9.3*, compared to the true angle of 0*. 

All the tests gave equally good predictions as the examples discussed above. Several 
conclusions can be drawn from the tests. First, there is no advantage and there is possibly a 
disadvantage in using an absolutely uniform array of probes. For example, if the sensor positions 
had been staggered vertically from one probe to another as shown in Figure 5.2-4b, the uncertainty 
in the position of the flat static position of the surface could have been reduced substantially. This 
improvement can also be obtained for static liquid surfaces that are curved, such as occur in low 
gravity. Second, the sensors in an array should be concentrated near the locations of the liquid 
surfaces used in the tests, rather than distributed uniformly along the entire depth of the tank. Third, 
the radial spacing of the probes should be no more than about 0.3R o - 0.4R o to obtain an adequate 
resolution of wave shape, and a closer spacing would be desirable near the wall to resolve the contact 
angle. Fourth, unless the direction of tank excitation can be fixed in advance and the array aligned 
with that direction, at least two probe arrays are required to resolve the line of peak slosh amplitude; 
these arrays cannot be at right angles. 
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5.3 Load Cell Requirements 

Measurement of liquid slosh forces is an essential requirement for both of the flight 
experiments. Therefore, a preliminary study was conducted to develop requirements for a force 
measurement system, and feasibility of a design for it. 

The peak slosh force at resonance can be estimated from: 

F - m,cfc 0 <2 (5.3-1) 

The slosh mass m, and natural frequency oo n for a spherical tank can be obtained from the results 
given in Section 4.0. The amplitude^ of the tank oscillatory motion and the resonance magnification 
factor Q can be estimated for linear sloshing to be: 

;t 0 = 0.05 R 0 (5.3-2a) 

Q = 25 (5.3-2b) 

For Bo ~ 0, the slosh force for the LH 2 spherical tank of the COLD-SAT flight experiment is thus 
estimated to be about 0.0006 lb (0.00025 N)). For the Shuttle flight experiment using water as the 
test liquid, the slosh force is about 0.0003 lb (0.0014 N) for the 13 in (33 cm) tank and 0.00018 lb 
(0.0008 N) for the 7 in (15 cm) tank. Compared to ordinary dynamics experiments, these are 
extremely small forces with extremely long (almost d.c.) time variations (see Tables 3.2-2 and 
3.3-2). 

A survey of commercially-available force transducers that may be suitable for flight 
experiments is summarized in [10]. The highest sensitivity listed in this survey is 500 mV/lb; thus, 
voltage variations of about 0.05 mV must be detected to measure the anticipated slosh forces. 
Custom-designed transducers using semiconductor strain gauges have therefore been used in 
previous laboratory studies of miniature slosh tanks [7, 27], and forces as small as 0.0001 lb [0.00044 
N] were measured reliably. Figure 5.3-1 shows a schematic of how these kinds of sensitive strain 
gauges can be used in a dynamometer to measure slosh forces, in this case in the plane of the figure; 
a similar set of dynamometers at right angles is required to measure slosh forces perpendicular to 
the plane of the figure. Note that the dynamometer is rigid in bending and actually senses the tensile 
and compressive strains in the legs. The data acquisition system used in [7, 27] cancelled the inertia 
force of the empty tank electronically, so the force that was detected was only the oscillating slosh 
force; this method of measurement increased the force sensitivity of the dynamometers by several 
orders of magnitude. Further optimization of the design would improve the force sensitivity even 
more. It is concluded, therefore, that a force measurement system for the flight experiments can 
be achieved. The specifications for such a system are listed in Table 5.3-1. 
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Both in-line and cross-axis response must be measured by means of an X-Y coordinate 
system. This is necessary in order to determine whether rotary slosh occurs. Static longitudinal 
dead load cancellation must be included in order to concentrate on the accuracy of the dynamic 
load measurement. Similarly, cancellation of empty tank dynamic load is desirable for the same 
reason. The requirement for the very low frequency range has already been mentioned, and will 
be aggravated by drift in the electronic circuits. Finally, static lockout will be required to avoid the 
effects of dead loads during high-g phases of the flights (both static and dynamic). Although all 
these requirements are extremely severe in terms of comparative requirements for high-g testing, 
it was concluded that a design which uses existing state-of-the-art technology is entirely feasible. 

TABLE 5.3-1 DYNAMOMETER REQUIREMENTS 

• lCF* to 10' 3 lbs. Range 

• In-Line and Cross-Axis Response 

• Static Longitudinal Dead Load Cancellation 

• Empty Tank Dynamic Load Cancellation 

• DC to 1 Hz Frequency Range 

• Static Lockout 
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6.0 CONCLUSIONS 
6.1. Flight Experiments 

This study of liquid slosh dynamics and control has shown that there is a real need for flight 
experiments to gain fundamental understanding about the phenomena that dominate liquid motions 
in low gravity and to improve and validate analytical models. Many planned space missions depend 
upon the ability to predict and control liquid motions, using information that can only be gained 
from realistic flight experiments. To that end, two flight experiments were defined: an experiment 
using liquid hydrogen in the spherical and cylindrical "receiver" tanks of NASA’s proposed 
COLD-SAT test bed satellite; and an experiment designed for the Shuttle middeck using water or 
other nonhazardous liquids. The types of liquid motions to be investigated include small-amplitude 
sloshing, nonlinear or rotary sloshing, and the decay of large-amplitude motions to long-lived, 
small-amplitude sloshing. To produce these motions, the test tanks will be subjected both to transient 
impulsive accelerations and to oscillatory accelerations sustained for up to eight to ten cycles. Both 
experiments will investigate the range of Bond numbers near zero. For the COLD-SAT experiment, 
the Bond number for each tank can be varied over about a factor of four by employing the satellite 
thrusters, while for the mid-deck experiment, the Bond number is constant. The specifications and 
data requirements for the flight experiments are within the capabilities of each respective carrier. 

Analytical efforts and ground-based experiments were conducted in support of the flight 
definitions. The conclusions from these efforts are summarized separately below. 

6.2 Analytical Studies 

6.2.1 Small Amplitude Linear Sloshing 

The computational fluid dynamics codes FLOW-3D and NASA-VOF3D were used to 
simulate low gravity sloshing in spherical tanks. Neither code was able, however, to make a realistic 
simulation, primarily because of their difficulty in modeling the "free" contact line condition (i.e., 
constant contact angle) and in accurately representing surface tension forces. Therefore, a new 
analysis of linear sloshing in low gravity was developed in an integral-minimization form; this 
formulation permits the relatively simple surface physics assumptions used in the analysis to be 
modified readily when better knowledge becomes available from flight experiments, which is 
thought to be a distinct advantage. Several errors in previous analyses of the slosh force were also 
corrected. "Trial" solution functions were computed by an innovative use of a finite-element 
structural code, which were then used in the integral-minimization technique to compute the sloshing 
frequencies, forces, and mechanical model. Since the method is based on a finite element structural 
code, it can be adapted readily by other investigators. 
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Examples were computed for a spherical tank to demonstrate the use of the methods. The 
examples employed fill levels of 25%, 50%, and 78%, for Bond numbers of 1 and 2. Methods were 
developed to "start" the finite element structural simulation such that rapid convergence of the 
numerics is ensured. The results of the examples show that the non-dimensional slosh frequency 
decreases as the Bond number decreases, in agreement with the limited experimental data available, 
although the frequencies were consistently slightly larger than previous finite difference simulations. 
The slosh masses of the mechanical model were predicted to be substantially smaller than the 
corresponding high-g slosh masses. 

6.2.2 Rotary Slosh Model 

An existing pendulum analogy of rotary sloshing was investigated more thoroughly and 
extended in a preliminary way to low gravity conditions, in order to estimate the importance of 
nonlinear effects in low gravity sloshing. 

A numerical scheme for integrating the equations of motion of the compound spherical 
and linear pendulum was developed. The scheme can be used effectively to establish all the pertinent 
model parameters for rotary slosh and is sufficiently general to be applied to any slosh test in which 
the cross-axis and in-line slosh forces (i.e., effective weights of the pendulum masses) can be 
measured. Some initial guessing of damping parameters is required, with help from the harmonic 
balance approximate model, to start the process. The guessing process also is helped by the fact 
that, at a given frequency point, an increase of both the circumferential damping and the radial 
damping tends to increase the cross-axis weight produced. Since damping increases in low gravity, 
this result indicates an increased tendency for cross-axis forces in low gravity, which thus impacts 
the design of control systems. 

The rotary slosh problem is complex in the sense that the mechanical model parameters 
vary with frequency throughout the pertinent range. Hence, no one set of parameters can be claimed 
to represent the compound system. Because of this complexity, the questions remains of what is 
the best approach to studying the system response to transient inputs. If one set of parameters must 
be selected, those associated with a > 1.0 are the most appropriate. This is the frequency range in 
which the maximum rotary motion occurs and the largest rotary forces. The harmonic balance 
approximation provides good estimates of the system parameters when the response is nearly 
circular. The next logical step is to make a numerical study of the transient response of the compound 
pendulum model for low gravity conditions. It is also appropriate to determine whether the addition 
of other mechanical components to the model might eliminate some of the coupling problems as 
well as the variation of the system parameters with frequency. Finally, although only 
counter-clockwise rotational motions were used in the present modeling study, clockwise solutions 
alsoexist for each case. The direction of rotation that actually occurs depends on the initial conditions 
and any lack of true rotational symmetry in the physical system. 
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6.3 Ground-Based Experimentation 

6.3.1 Liquid-Vapor Interface Sensors 

Miniature electrical resistor sensors and, possibly, fiber optic sensors can be used to 
determine the location of a discrete point on a liquid interface. An array of such sensors can therefore 
be used in flight experiments to track and "visualize" the motions of the liquid surface. The frequency 
response of these sensors is substantially better than the relatively slow frequencies of the liquid 
surface motions expected in low gravity. Ground testing, using non-flight sensors, was employed 
to determine the spatial distribution of such arrays of sensors that are required to resolve low gravity 
surface configurations and slosh wave shapes. Since the sensors give only a "wet to dry" or "dry 
to wet" indication, an interpolation scheme was also developed as part of the testing to convert this 
digital data from the sensors to a continuous time history of the motion of each discrete surface 
point. The test results proved that the free surface configuration and the slosh wave shape could 
be resolved satisfactorily by arrays that were spaced radially at distances of no more than 15% to 
30% of the tank radius, and spaced vertically at about 10% of the tank radius. Better resolution 
could be obtained with denser arrays, especially near the tank walls. In any case, the number and 
spacing of the sensor arrays can be optimized for specific liquid levels. It was concluded that ( a ) 
the frequency response of available sensors is more than adequate for flight experiments, and (6) 
the required number and distribution of the sensors does not appear to be impractical. 

6.3.2 Slosh Force Transducers 

The slosh forces that must be measured are in the range of 0.0006 lb (0.00025 N) for the 
COLD-SAT flight experiment and 0.0002 lb (0.0009 N) for the middeck flight experiment. The 
natural period of the sloshing is about 800 sec for the COLD-SAT experiment and 100 sec for the 
middeck experiment. Such small forces with such relatively long periods cannot be measured with 
conventional laboratory transducers. Nonetheless, a survey of commercially-available transducers 
and of dynamometers used previously in laboratory testing of small Bond number sloshing did 
reveal that the force measurement requirements can be satisfied. The most promising of the available 
methods used semiconductor strain gages in a tension-compression arrangement, coupled with 
electronic cancellation of non-slosh forces, to measure slosh forces as small as 0.0001 lb 
(0.00044 N). Further optimization of this design, which is inherently rugged and capable of 
supporting large dead loads, should satisfy the requirements of the flight experiments. 
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Figure 3.3-1 Conceptual Design of Shuttle Flight Experiment Package 
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Figure 3.3-2 Slosh Characteristics of a Half-Full Spherical Tank 
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Figure 3.3-3 Block Diagram of Shuttle Flight Experiment Instrumentation and Hardware 
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(b) Contact line projections on the tank wall 


Figure 4.4-2 Details of Low-G Slosh Analysis 
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Figure 4.4-3 
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Figure 4.4-6 Slosh Wave Shape for 50% Full, Bo = 


s 


Figure 4.5-1 
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Figure 4.5-3 Steady Deflection Angie 
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Figure 4.5-4 Position Lag Angie 
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Figure 4.5-5 Spherical Pendulum Mass Fraction 
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a. Space Plane (£0, <f>) 



b. Time-Phase Plane 


Figure 4.5-6 Relationship Between Physical Space Plane and Time-Phase Plane 
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Figure 4.5-7 Polar Plot / sin 6 - Inches, cj> - Radians 
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Figure 4.5-8 Deflection Angle, 0 - Radians 
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Figure 4.5-9 Circumferential Angle, <|> - Radians 
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Figure 4.5-10 Position Angle, <j> 0 , - Radians 
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Figure 4.5-13 Polar Plot (a = 1.000) 



Figure 4.5-14 Polar Plot (a = 1.107) 



Figure 4.5-15 Influence of Radial Damping on Character of Response 
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Figure 4.5-17 Spherical Pendulum Deflection Corresponding to Rotary Slosh Amplitude 









Figure 4.5-18 Dependence of Spherical Pendulum Response on Frequency and Damping 
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Figure 5.2-1 Schematic of Fiber Optic Liquid-Vapor Sensor 



Figure 5.2-2 Slosh Test Apparatus for Liquid- Vapor Interface Sensors 
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(b) Staggered array 


Figure 5.2-4 Probe Arrays for Liquid- Vapor Interface Sensors 
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(a) Probe R3 
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(b) Probe R2 

Figure 5.2-5 Interpolated Test Data and Predicted Wave Shape for Four-Sensor Probe Array 
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(d) Wave shape 


Figure 5.2-5 (Concluded) 
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Figure 5.2-6 Interpolated Test Data and Predicted Wave Shape for Six-Sensor Probe Array 
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Figure 5.2-6 (Concluded) 
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Figure 5.3-1 Strain-Gage Slosh Force Dynamometer 
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1.0 INTRODUCTION 


This document describes the requirements for the Liquid Dynamics and Slosh Control 
Experiment, a Class I experiment, to be performed on COLD-SAT. The purpose of 
the experiment is, primarily, to gain a physical understanding of the motions 
of liquids having a free surface in a low gravity environment, including the 
acquisition of detailed data for several representative cases; this understanding 
can then be used to develop and verify Improved analytical models. A secondary 
purpose is to characterize the dynamic interaction of the liquid motions with 
a representative spacecraft (COLD-SAT). Included in the document are some 
background information on low-g liquid dynamics, a statement of the objectives 
of the experiment, a description of the physical parameters and processes to be 
investigated, a description of the experiment, a preliminary test matrix, the 
experimental procedures, and the data required from the experiment. 

Because free-surface motions in a tank depend strongly on container geometry, 
liquid fill level, liquid physical properties, ambient gravity level , and container 
motion, it is not practical to determine the dynamics of the motions by low-gravity 
experimentation for every mission anticipated by NASA and DOD. Instead, analytical 
and numerical models must be used for most of the missions. But these models 
must be validated by comparison to a reference set of actual low-gravity data 
for a few representative cases. The scope of the COLD-SAT experiment has been 
formulated with this goal in mind. 

2.0 BACKGROUND 

Liquid dynamics in the tanks of space vehicles has long been recognized as being 
important to stability and structural loading. NASA's plans for an ambitious 
Space Exploration Initiative to the Moon and Mars Involve vehicles that will 
transport and store enormous quantities of cryogenic propellants in space. 
Additionally, some kinds of space-based optical systems (e.g., the Strategic 
Defense Initiative) also Involve large mass fractions of liquids. Tank sizes 
for these and other space vehicles range from a meter or so In diameter to several 
tens of meters. The control, pointing, and docking of spacecraft containing 
such large masses of liquid, as well the transfer of liquid between spacecraft, 
is critically dependent upon understanding and controlling the motions of the 
liquids. 

The dynamics of liquids In tanks in normal gravity is well understood, and 
analytical, numerical, and scale-model test methods have been established to 
treat these "high-g" problems [1, 2]. However, low-gravity free-surface liquid 
motions are not nearly so well understood. The motions are dominated by surface 
physics effects that cannot be Investigated realistically by ground testing in 
normal gravity. Some Information is available from drop tower and zero-g tests 
[e.g., 3] but all such studies to date have employed small tanks and non-cryogenic 
liquids. It appears from this limited amount of data that, for reasons that are 
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not well understood, slosh damping is larger in low-gravity [3, 4] and nonlinear 
effects such as rotary sloshing are more prominent [5, 6] . As additional evidence 
of the importance of surface physics, ground-based slosh tests that simulate 
low-gravity indicate that the response of the liquid depends strongly on the 
motion of the liquid-tank contact line [e.g., 7]; for example, the natural 
frequency and damping varies by more than a factor of two between a "free" contact 
line and a "stuck" contact line condition. Consequently, this COLD-SAT experiment 
is needed to study surface motion of cryogens in tanks of moderate size, in an 
actual low-gravity environment for a sufficiently long test duration. 

3.0 TECHNOLOGICAL OBJECTIVES 

The specific objectives of this technology experiment are: (1) determine and 
understand the liquid motion resulting from typical maneuvers of spacecraft; and 
(2) characterize the Interaction of the motion with spacecraft for a typical 
spacecraft (COLD-SAT). These objectives will (a) provide physical understanding 
about the motions of cryogens In low-gravity; (b) supply data from which to 
establish the adequacy of existing analytical and numerical models of such 
motions, and indicate where Improvement Is needed; and (c) yield typical data 
on the damping of low-gravity motions from which analytical and empirical 
correlations can be developed. 

The physical processes to be investigated include: 

Static liould configuration - The configuration of the liquid in the test 
tanks will be monitored under ambient on-orbit conditions. 

Liquid response is various discrete accelerations - The liquid will be 
oriented to give a specific Initial orientation. Impulsive and periodic 
accelerations of selected amplitude, frequency, and duration will be 
applied and the free surface response monitored. 

Fffpctiveness of slosh baffles - The cylindrical receiver tank will contain 
a single ring baffle to demonstrate the damping of liquid motions in low 
gravity. 

4.0 JUSTIFICATION 

The motion of contained liquids has a profound Influence on the dynamics and 
control of space vehicles, the transfer of liquids between vehicles, and the 
docking of one vehicle to another. Future space missions will carry much larger 
quantities of liquids, primarily cryogens, than is common now, and the mission 
performance requirements will be much more demanding. As an example of the 
problems that must be solved, space-based telescopes and strategic defense 
satellites must be pointed to an angular accuracy of the order of 0.0005* and 
that accuracy must be maintained even during tracking maneuvers; since the 
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tracking maneuver will set into motion the contained liquids, one can imagine 
that the ability to understand, predict, and control the dynamics of the moving 
liquid will be critical. 

Although some spacecraft maneuvers will certainly create large motions of the 
contained liquid, this experiment will concentrate on liquid motions that are 
localized about the initial position of the liquid; such motions are usually 
called "sloshing." There are several reasons for limiting the experiment in 
this way. First, large motions are dominated by liquid inertia and so can be 
modeled more easily than motions dominated by little-understood surface physics. 
Second, large motions eventually decay to smaller amplitude sloshing motions. 
Third, slosh motions can be excited by a variety of typical control maneuvers 
of spacecraft and are thus important in their own right. 

5.0 ANALYTICAL MODELS 

Low gravity fluid dynamics in "bare" tanks are currently modeled either by general 
purpose computational fluid dynamics (CFD) codes, such as NASA-V0F3D or FL0W-3D 
[8, 9], or by special-purpose linearized analyses [10], At present, however, 
only the CFD codes can model slosh baffles and other internal hardware. Both 
methods employ assumptions about the surface physics that either need to be 
validated or improved by in-space experimentation with cryogens. 

The parameters needed for the models include tank shape, liquid fill level, and: 

Bond number Bo=g eJf a 2 /$ 

Static contact angle 0, 

where g eff is the effective "gravity" or linear acceleration of the spacecraft, 
a is a representative dimension of the tank, and p - a/p is the specific surface 
tension of the liquid. The tank shape and the direction of g eff with respect to 
the tank axis must also be specified, and, for Investigations of nonlinear 
effects, the tank motion as well. 

Bond numbers Bo < 10 or so represent "low" gravity conditions, while Bo « 1 
represent "micro" gravity conditions; Bo < 1 is representative of a large tank 
(> 1-m in diameter) in orbit. Whenever the contact angle of the liquid with the 
tank is small (as it is for cryogens and tank materials of interest), the free 
surface of the liquid is highly curved in low gravity; the extreme case of Bo 
- 0 and a zero-degree contact angle results in a free surface that is spherical. 
The stability of the free surface is also a function of Bo (based on the disturbance 
acceleration) and contact angle. 
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The models also require a relationship that specifies how the contact angle 
changes as a function of the contact line velocity when the liquid surface is 
in motion. Currently, the only relation that can be employed in the models is 
that the contact angle remains constant at its equilibrium value, regardless of 
contact line velocity. Although this assumption appears to agree with most (but 
not all) of the available small-tank, drop-tower test results, fundamental studies 
of the spreading of liquids on surfaces indicate that it may be an over-idealization 
[11]. This COLD-SAT experiment is designed to establish better contact angle 
- contact line velocity relations for cryogenic liquid and tank materials of 
Importance. If the results of the experiment show that an assumption of constant 
contact angle is not valid, the models must be modified. Other surface physics 
phenomena may be important and therefore will need to be included in the models 
(e.g., surface tension hysteresis) but the nearly complete lack of low-gravity 
slosh data has so far prevented the establishment of all potentially important 
parameters . 

There are currently no models of the damping provided by slosh baffles in 
low-gravity. The models that will be evaluated are modifications of high-g 
results [1]. The models of low gravity viscous, or "bare" tank, damping have 
not been well validated; for example, some correlations developed from drop-tower 
and simulation tests show that the viscous damping ratio y is a function of Bond 
number as well as viscosity: 

y= A (v/f„a 2 ) 1/2 [ 1 + C (Bo )""] 

where v is the kinematic viscosity and f„ Is the slosh natural frequency. These 
correlations were developed from tests with Bo > 1 and may overestimate the 
damping significantly when Bo < 1 (because of the negative exponent on the Bo 
term). Improved correlations will be developed from the results of this COLD-SAT 
experiment. 

6.0 EXPERIMENT REQUIREMENTS 

6.1 Description of Experiment 

To investigate the Influence of tank shape, slosh experiments will be conducted 
in both receiver tanks. The cylindrical receiver tank will be fitted with a 
single ring-baffle near Its midpoint, but the spherical receiver tank will not 
contain any specific anti -slosh devices. The tank support structures will 
incorporate load cells to monitor the forces and torques exerted on the tanks 
by the sloshing liquid. Liquid-vapor sensors in the tanks will be used to monitor 
the static liquid configuration and the motion of the free surface. The COLD-SAT 
propulsion system will be used to provide specific linear accelerations ( g eff ) 
and discrete disturbance accelerations. 
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These COLD-SAT experiments can be performed as opportunities (which are plentiful) 
when the liquid fill volumes of interest are available during other tests. Some 
of the specified disturbance accelerations may also occur during other testing. 

6.2 Key Parameters 

The key parameters expected to influence the low-g slosh experiments are (1) 
liquid surface tension, density, and viscosity (which are functions primarily 
of liquid temperature for a given liquid), (2) tank shape and internal hardware, 
(3) fill level, and (4) acceleration environment. It is not proposed to vary 
the liquid temperature systematically, so the effect of liquid properties will 
not be investigated except as the temperature may change for other reasons from 
test to test. The influence of tank shape will be investigated by tests with 
both receiver tanks, one of which is basically a cylinder and the other a sphere. 
The only internal hardware that will be specifically investigated will be a slosh 
baffle (an annular ring) in the cylindrical tank to investigate the performance 
of a typical baffle in low gravity when the static liquid surface is near the 
baffle; other static liquid levels will be chosen to simulate “bare wall" 
conditions in the same tank. The influence of liquid level will be investigated 
by conducting tests, to the extent possible, over a variety of fill levels. The 
influence of effective gravity will be investigated by conducting tests at several 
thrusting levels, as described later in Section 6.5; the acceleration will be 
directed along the symmetry axis of the tanks to facilitate comparison with 
analytical /numerical models. Perturbation accelerations to excite sloshing will 
include impulsive and periodic motions of the spacecraft. The duration of the 
periodic acceleration will be varied, depending upon the objective of the test, 
as described later in Section 6.5. 

6.3 Measurements 

The static and dynamic configuration of the free surface is important in 
understanding and correlating the experimental results; for example, the slosh 
wave shape near the wall will indicate whether the contact angle remains constant 
during sloshing or, alternatively, that the contact line is "free." Liquid-vapor 
sensors mounted in the receiver tanks will therefore be used to monitor both the 
static and dynamic configurations and the frequency of the surface motions. To 
accomplish this, the sensors must be arrayed vertically and radially, as 
illustrated in Figure 1; ideally, the sensors should be contained in two orthogonal 
arrays. Ground-based experiments indicate that the free surface shape can be 
resolved accurately when the radial spacing of the sensors is as large as 20 % 
of the diameter (i.e., spaced at 10%, 30%, 50%, 70%, and 90% of the diameter) 
and the vertical spacing is as large as 10% of the diameter. More widely spaced 
intervals can be used without sacrificing accuracy in the frequency estimation 
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Figure 1. LI 


(which Is determined by the data sampling rate) but the surface shape estimation 
will be degraded. Note that the vertical spacing of the sensors can be concentrated 
at the static liquid levels employed in the tests, and, if necessary, the radial 
array locations can be reduced to 10%, 30%, and 50% of the tank diameter. Because 
of the low frequencies of the anticipated slosh motions, the required frequency 
response of the liquid-vapor sensors is essentially d.c. Of more importance is 
the time required for a sensor signal to change from “dry" to "wet" and from 
"wet" to "dry" as the liquid surface passes the sensor location; ground-based 
experiments Indicate, however, that the signal changes abruptly over a time 
period that Is very small compared to the time between successive passages of 
the liquid surface; thus, the sensors should function more than adequately. 

Knowledge of the overall slosh wave shape and frequency, as acquired by the 
liquid-vapor sensor array. Is sufficient to examine the predictive accuracy of 
analytical models. More fundamental information may be needed to identify weak 
assumptions in the models. For example, the relation between the dynamic contact 
angle at the wall and the contact line velocity can potentially have a significant 
influence on the slosh frequency and force, even though existing models assume 
that the dynamic contact angle remains equal to the static contact angle. The 
liquid-vapor sensor array Illustrated In Figure A.l will not be adequate to 
acquire detailed data about the wave shape near the wall. To resolve the wave 
shape near the wall in detail, and thus to infer the behavior of the contact 
angle as a function of time, will require a denser radial and vertical array of 
sensors near the wall for at least one selected liquid level. These sensors 
should be spaced symmetrically at about 0.02R, above and below the free surface 
over a total vertical distance of about O.lRo. and at least two such arrays should 
be spaced radially at Intervals of 0.02R, from the wall. Because of the "bent 
over" geometry of the tank-liquid intersection for spherical tanks, the dense 
arrays are most readily made applicable to cylindrical (straight wall) tanks, 
although a dense array of sensors positioned along a radial line could be used 
for spherical tanks. In addition, the array installation must not significantly 
interfere with the slosh wave motion. It is understood that the use of such 
dense arrays will require substantial data acquisition rates and may interfere 
with other experiments. 

The support structure of the receiver tanks must be designed to accommodate load 
cells with which to measure slosh forces. The slosh force measurements can also 
be used as a means of Independently determining slosh natural frequencies and 
damping, in the event of difficulties in interpreting the liquid-vapor sensor 
data for a particular test. The magnitude of the slosh forces that must be 
measured is a function of Bo and the amplitude of the Imposed disturbance 
acceleration. For the COLD-SAT receiver tanks and the proposed excitation 
amplitudes, the expected force amplitudes range from about 0.0001 lbs to about 
0.001 lbs. The required frequency response of the load cells also depends on 
Bo but is essentially d.c. Such small forces can be measured reliably only by 
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eliminating (or otherwise compensating) the large non-sloshing dead-load signal 
from the dynamic slosh force signal. This compensation can be accomplished 
electronically before the signal is transmitted to the data acquisition system, 
by methods used previously In ground tests [e.g., 4], 

The liquid temperature near the free surface and tank pressure must be measured 
at the initiation and conclusion of each test, to determine the liquid properties. 
The accuracy of the required pressure measurement is ±5 psla, and the accuracy 

of the required temperature measurement is ±1°R. The temperature determination 
must be made at enough locations to estimate a representative average temperature. 

TABLE 1 

MEASUREMENT REQUIREMENTS FOR LIQUID DYNAMICS EXPERIMENT 


Parameter 

Range 

Accuracy 

Instrument 

Liquid temperature 

20 - 50’ R 

± l’R 

COLD-SAT temp, 
sensors 

Tank pressure 

5-50 psla 

± 5 psla 

COLD-SAT pressure 
sensors 

Liquid interface 

Cylindrical tank: 

± 5% axial 

Liquid/vapor sensors 

position 

0 - 4 ft axial; 

0 - 1.3 ft. radial 
Spherical tank: 

0 - 1.3 ft 

± 10% radial 

with 2 - 5 second 
response 

Slosh force 

0 - 0.0001 lb 

± 2% 

Load cells on tank 
support structures 

Slosh frequency 

0.001 - 0.01 hz 

± 1%. 

Load cells and 
liquid/vapor sensors 

Steady linear 
acceleration 

8p.g - lOOpg 

± 1% 

COLD-SAT acceler- 
ometers or compute 
from thruster firing 
histories 

Perturbation 

acceleration 

Impulse: A t m 8pg 
Periodic: A p -8pg 
V 75 - 125 sec 

± 1% 

COLD-SAT gyros and 
accelerometers, or 
compute from thruster 
firing histories 


The effective gravity acceleration and the magnitude and history of the disturbance 
accelerations will be measured by the satellite accelerometers and gyros. The 
effective gravity amplitude should be measured to ±1%. Ideally, the perturbation 
accelerations should also be measured to within ±1% In amplitude and continuously 
in time. If these requirements cannot be met, the steady and perturbation 
accelerations can be computed from the history of the specified thruster firings. 

Table 1 summarizes the measurement requirements. 
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6.4 Hardware Requirements 


The special hardware items required for the experiment are: (1) annular ring 
baffle mounted in the cylindrical receiver tank; (2) load cells mounted on the 
support structures of both receiver tanks, and (3) liquid-vapor sensor arrays 
in both receiver tanks. 

The recommended width of the ring baffle is 10% of the tank radius. The baffle 
should be located at an axial position such that (a) the liquid free surface for 
40% filling is slightly above the baffle and (b) no part of the free surface is 
pierced by the baffle. This location will maximize the damping at the 40% filling 
level while permitting higher filling levels to be investigated without significant 
damping from the baffle. 

The load cells must be mounted on the tank support structures at locations where 
the lateral forces can be sensed. 

The liquid- vapor sensors must be arrayed In sufficient numbers to determine the 
liquid interface position with the accuracy indicated in Table 1. 

6.5 Procedures 

The liquid dynamics experiments will be performed primarily when opportunities 
arise. Whenever the satellite is maneuvered for other experiments or for 
operational reasons, the data necessary to define the quantity of liquid in each 
receiver tank, the disturbance, and the liquid-vapor sensor responses can be 
acquired. Nonetheless, several tests using periodic excitation probably will 
not occur as opportunities and must be conducted specifically. 

The tests identified in Table 2 will be performed when the desired liquid levels 
(± 5%) are reached over the course of testing. During the tests, the satellite 
attitude control thrusters will be operated with specific commands to obtain the 
desired perturbation accelerations. The initial orientation of the liquid will 
be established by applying a settling acceleration with the thrusters, and this 
level of steady acceleration must be maintained for the duration of each test. 
The motion of the liquid surface will be monitored by the liquid-vapor sensors 
and the slosh forces by the load cells throughout each test. 

The steady accelerations listed in Table 2 were selected to conform to the nominal 
levels that can be obtained by firing a single engine or a combination of engines; 
the exact acceleration level Is not Important so long as it Is recorded for later 
data analysis. The perturbation accelerations are nominal also, and correspond 
to the firing of the appropriate attitude control thrusters or, depending on the 
COLD-SAT design, a gimbaled engine; again, the exact levels of the accelerations 
are not important with the exception that they should be a small fraction of the 
steady acceleration. Impulsive perturbations are obtained by firing the engines 
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TABLE 2 PRELIMINARY TEST MATRIX 
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for a short time, and periodic perturbations by on-off firings for scheduled 
periods. The frequency of the periodic accelerations should be maintained as 
indicated in the table, since the intent is to excite near-resonance sloshing. 

The tests using impulsive excitation are designed primarily to determine the 
slosh natural frequency and damping. The tests using periodic excitation are 
designed either to determine slosh forces or to evaluate nonlinear effects, 
particularly the tendency for rotary slosh. 

6.6 Data Analysis 

Liquid temperature and pressure measurements will be used to determine the liquid 
surface tension, viscosity, and density from tabulated data. Liquid volume and 
liquid-sensor measurements before the perturbation acceleration is applied will 
be used to estimate the initial configuration of the liquid surface and the 
static contact angle. Accelerometer and gyro data acquired during the test (or 
computations based on the thruster firing histories) will be used to determine 
the imposed motion of the test tanks. 

For impulsive perturbations, liquid-vapor sensor measurements during the sloshing 
will be analyzed to compute (1) slosh wave shape and amplitude as a function of 
time, (2) slosh natural frequency, and (3) slosh damping (from the decay of the 
slosh wave amplitude). If a dense array of sensors is used near the wall as 
discussed in Section 6.3, contact angle and contact line velocity will also be 
computed as a function of time; if not, these quantities will still be estimated 
but the resolution is not expected to be sufficient for fundamental studies. 
The load cell force histories will be analyzed to confirm the slosh natural 
frequency and damping data and to compute the slosh force. For periodic 
perturbations, the data analysis will be similar to the impulsive acceleration 
tests, with the exception that neither the liquid-vapor sensor nor the load cell 
measurements will provide damping data. The load cell data, in conjunction with 
damping data from the impulsive tests for the same Bond number and liquid level, 
will be analyzed to determine the effective mass of liquid participating in the 
sloshing. 

Eventually, all the analyzed test data will be used to compare with predictions 
from the analytical /numerical models for the same Bond number, fill level, contact 
angle, and excitation. 
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Table 3 NOMENCLATURE 


a tank radius, cm 

A t amplitude of impulsive acceleration, g's 

A p amplitude of periodic acceleration, g's 

Bo Bond number, g^a 2 / P 

/„ slosh natural frequency, hz 

g eff effective gravity or steady linear acceleration, cm/sec 2 

p specific surface tension, a/p, cm 3 /sec 2 

v kinematic viscosity, cmVsec 

0, static contact angle, degrees 

p density, g/cm 3 

a surface tension, dynes/cm 

Tj length of impulsive acceleration, sec 

Tp period of periodic excitation, sec 

x, slosh period, 1 //„, sec 
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APPENDIX B 
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B.l Solution for Equilibrium Surface Shape 

The free surface of a fluid in a tank is defined by the following differential equation in the 
normalized surface coordinates introduced in Section 4.0. 


c£ZdR 
ds 2 ds 


c£R_dZ 
ds 2 ds 


B °(z -z 0 ) + If - S. = 


Five boundary conditions are required to solve this equation, 
there are four conditions which are specified, 


0 (B.l-1) 

At the centerline of the tank 


R= 0, z = -z 0) f = 0, ^= 1 , 

ds ds 


at s =0 


(B.l -2) 


The fifth condition is defined by the contact angle between the fluid and the tank wall, 


dZ o 
— = cos 0. 
ds 


at s = s woU for cylindrical tank 


(B. 1-3) 


: 0, + sin 1 

'ZoX 

c 

J. 


at s = s waU for spherical tank 


Because Equation (B.l-1) is highly nonlinear, a computational approach to its solution is 
used. This equation can be reduced to a set of coupled first order equations with the following 
approach. Consider the relation between the surface coordinates. 


dRY 

ds j 


+ 


f dZ Y 
j 


(B.l-4) 


Solving this relation for dR/ds and differentiating. 


dR = \ JdZYY 

ds L t ds J. 


d 2 R d 2 ZdZ(dR_Y 

ds J 


ds 2 ds 2 ds 

Substitute this into Eq. (B.l-1) and rearrange to yield, 


d^Z 

ds 2 


dR 

ds 


\+Bo{Z-Z 0 )-~ 
R ds 


(B.l -5) 
(B.l-6) 


(B.l-7) 
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Now solve Eq. (B.l-4) for dZ/ds and differentiate, 

r dR 


dZ_ 

ds 

d*Z 

ds 2 


1- 
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\2~l 2 


ds 


d 2 R dR 


ds 2 ds 


dZ 

\ ds J 


V* 


Again, substitute this relation into Eq. (B.l-1), 


S = -R + B0(Z-4)]f ^ 


Four first order differential equations can now be formed. 


(B.l-8) 

(B.l-9) 
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ds 
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dZ^ 

ds 


dR x 

ds 


Xg + Bo (Z x Z 0 ) 


J_dZi~ 
R x ds 


dR x 

ds 


= R 
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(B.l-12) 

(B.l-13) 


dRr 

ds 
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These four equations can be solved with the use of a standard fourth-order Runge-Kutta 
integration algorithm. The four initial conditions at s=0 are specified to begin the solution. The 
solution is obtained with various values of until the contact angle condition is satisfied. This 
method is similar to the "shooting" technique used to solve the Blasius equation for a laminar 
flat-plate boundary layer profile. 

It should be noted that Eqs. (B.l-12) and (B.l-14) contain a singularity at R=0 which must 
be resolved before the solution can proceed. First, take the limit of Eq. (B.l-12), 


{ 


lim 

t -»o 


V 


dZj 

ds 


A 

) 


\ + lim 


J_dZ, 
R x ds 


\ 

J 


(B.l-15) 


Using L’Hospital’s rule. 
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(B. 1-16) 


or. 




lim 
/ — > 0 


dZi 

~ds 


= K - 


lim 

s ->0\ 1 is* J 


lim 

x ->0! 


dZ, 

~ds 


- 5 ^ 


Similarly, Eq. (B.l-6) can be used, 

d 2 Z,dZ x (dR^ 


lim 

x — »0l 


r dR A 


ds 


4 -^ 

x-» 0 ^ ds J 


= limi 
x-»o[ ds ds 

= 0 


ds 


x-fo^ ds 



(B.l-17) 


(B.l-18) 

(B.l-19) 


A computer program was developed for solving the system of equations, (B.l-1 1) - (B.l-14). 
A listing of this routine follows. This routine requires the user to manually iterate on the value of 
Xo until the desired solution is achieved. This particular version of the code will execute on an 
IBM-PC when compiled with Microsoft FORTRAN 5.0. Very minor modifications are required 
to execute the code on other platforms. 

An automatically iterating version of this program which uses a Newton-Raphson approach 
to converge on a value of X<, was also prepared. This approach is highly sensitive to the specific 
conditions and requires certain program modifications to "tune” the convergence rate for each 
application. So, it is not presented here. 

A sample input File for this program follows the program listing. This file is read with 
FORTRAN list-directed READ statements; so, the format is rather flexible. A descriptor of each 
value in the file can be included because of this choice of file reading technique. The values listed 
here should be contained in a file named GENSURF.INP. These values are appropriate for a 
spherical tank that is 75% filled for the case of Bo-1. A value of X 0 = 2.98068 (supplied by the 
user from the keyboard) will provides a contact angle of 5.0467* , which is within 1% of the desired 
values of 5°. 
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B.2 Solution for Slosh Potential and Mechanical Model Parameters 

The non-dimensional slosh velocity potential is defined as: 


& = I W (B.2-1) 

where the Ui are the trial solutions from the structural finite element simulation. The b { are the 
modal participation factors, or the eigenvectors, determined from the eigenvalue matrix equation: 

[A] {17} = Q 2 [fi]{17} (B.2-2) 

Equation (B.2-2) is derived from numerically integrating the minimum integral, Eq. (4.4-29), with 
the assumption that dUJdN = Q 2 (l +Bo)MU i on the free surface. The matrix A is: 


f f - 

“ 

r dR y 

Wmu, 

BoR 


J i 
0 1 


ydS j 


d 

BS 




-R 


(l_ J_ 

2 



(B. 2-3) 


and the matrix B is: 

jMUj/jdS 

0 


(B.2-4) 


A QUICKBASIC computer code "LOW-G.BAS" was written to integrate Eqs. (B.2-3) and (B.2-4) 
numerically. The code computes the required S-derivatives of Uj of [A] numerically by fitting a 
quadratic through the value of the central point and the two adjoining points; the points at the S = 
0 and S = S c are handled by special formulas. 

The code also solves the matrix equation (B.2-2) to determine the eigenvalue Cl 2 for the 
fundamental slosh mode and the corresponding b t eigenvectors. The slosh wave amplitude H c at 
the wall is then computed: 


H„ 


(l+Bo)M(S=S c ) 

Q 2 




(8.2-5) 


The slosh mass is computed by numerically integrating Eq. (4.4-28): 


npRo 



[(2K^? e - cos y.) cot y c ]H c 
( l+Bo)Q 


r, 



o 


(B.2-6) 


where the subscript w indicates that the parameter is evaluated at the tank wall. The normal derivative 
of O in the denominator of Eq. (B.2-6) is computed from analogy to Eq. (B.2-1) as: 


do 

dN 


MQ+Bo^bfoUi 


(B.2 -7) 


v m (V-| jarer* f \w . 'f **** 
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The other parameters of the mechanical model are computed from Eqs. (4.4-26a) and (4.4-25b). 
The QUICKBASIC program also computes these parameters. 

The code requires the following general input information: 

Bond number 
Filling percent 
Contact angle 

Exponent m and parameter e of Eq. (4.4-36) 

Number of nodes on the free surface along the line 0 = 0 
Number of trial modes (1, 2, or 3) 

The code also requires the following data for each of the finite element nodal points on the free 
surface along the line 0 = 0 

Nondimensional nodal coordinates R, S 

Nondimensional derivatives dR/dS, (fRIdS 2 , dZ/dS, and cfZIdS 2 

Trial potentials £/, for each mode 

and the data at each finite element nodal point on the tank wall along the line 0 = 0: 

Nondimensional nodal coordinates R, Z 
Trial potentials U ,• for each mode. 

The coodinate data and the various derivatives are available from either the finite element simulation 
model or from the code described in Section B.l. 

A sample set of input data is shown in the following pages for a 50% full tank and Bond 
number of one. As shown, the data is entered from the terminal. Any of the data for each screen 
can be corrected when the screen is displayed. The data is saved in two disk files COORD .DAT 
and MODE.DAT and the case can be run again, and the input data corrected if necessary, by 
indicating on the first screen that the data is to be read in from the disk files. The computed output 
is displayed on the screen and, if desired, printed. 

A listing of the code is also attached. The code is extensively commented so that the logic 
flow can be readily followed. 
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Hit RETURN after each entry 


Bond No 


. = 1 Filling % = |7oJ 

notion: Exponent = 2.5 


Contact angle at wall (deg) - 


Mass function: Exponent = 2.5 Constant - 1.1 

How many trial modes do you want to consider? 3 

How many nodal values are there at the surface and wall? 


Eigenvalue of trial modes = 0.456469 1.53750 2.88534 


DO YOU WANT TO INPUT NODAL DATA FROM THE (1) TERMINAL OR (2) FILES ? 


■p 








DATA FOR FREE SURFACE NODES 
Enter: dR/dS and dZ/dS derivatives for each node. 

U sa tha ENTER or DOWN ARROW koy to antar aach valua from tba £11*. 
Uaa tha UP ARROW to move back up through tha data. 



Noda # 

R' (S) 

Z' (S) 


Noda # 

R' (S) 

Z' (S> 

1 

79 

1 

0 

13 

1821 

.45919 

.88832 

2 

248 

.99662 

.08198 

14 

1900 

.33774 

.94123 

3 

416 

.9865 

.16369 

15 

1967 

.21826 

.97588 

4 

582 

.96951 

.24499 

16 

2069 

.10303 

.99465 

5 

745 

.94548 

.32566 

17 

2069 

-.00601 

.99996 

6 

904 

.91412 

.40545 

18 

2106 

-.10773 

.99415 

7 

1058 

.87508 

.48395 

19 

2135 

-.20132 

.97949 

8 

1206 

.828 

.5607 

20 

2157 

-.2865 

.95804 

9 

1347 

.77246 

.63503 

21 

2173 

-.36323 

.93166 

10 

1480 

.70807 

.70612 

22 

2184 

-.43184 

.9019 

11 

1604 

.63449 

.77293 

23 

2191 

-.4928 

.87012 

12 

1718 

.55154 

.83414 

24 

2195 

-.54656 

.83738 





25 

2197 

-.59386 

.80457 
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cm m 


DATA FOR FREE SURFACE NODES 
Enter: Potential values for mode number 1 

Us* the ENTER or DOWN ARROW key to enter each value from the file. 
Use the UP ARROW key to move back up through the data. 



Node f 

Potential 


Node # 

Potential 

1 

79 

0 

13 

1821 

1.0046 

2 

248 

.08704 

14 

1900 

1.1151 

3 

416 

.17291 

15 

1967 

1.2131 

4 

582 

.25714 

16 

2069 

1.3021 

5 

745 

.34031 

17 

2069 

1.3826 

6 

904 

.42302 

18 

2106 

1.4543 

7 

1058 

.50533 

19 

2135 

1.5171 

8 

1206 

.58777 

20 

2157 

1.571 

9 

1347 

.6702 

21 

2173 

1.6165 

10 

1480 

.75276 

22 

2184 

1.6537 

11 

1604 

.83602 

23 

2191 

1.6834 

12 

1718 

.92011 

24 

2195 

1.7069 




25 

2197 

1.7266 


Input Screen No. 5 


DATA FOR FREE SURFACE NODES 
Enter: Potential values for mode number 2 

Use the ENTER or DOWN ARROW key to enter each value from the file . 
Use the UP ARROW key to move back up through the data. 


xie * 

Potential 


Node # 

Potential 

79 

0 

13 

1821 

- . 0827 6 

248 

.17499 

14 

1900 

-.47634 

416 

.33763 

15 

1967 

-.89097 

582 

.47736 

16 

2069 

-1.3111 

745 

.58819 

17 

2069 

-1.7256 

904 

.6655 

18 

2106 

-2.1232 

1058 

.70452 

19 

2135 

-2.4927 

1206 

.70172 

20 

2157 

-2,8277 

1347 

.65315 

21 

2173 

-3.1221 

1480 

.55491 

22 

2184 

-3.3721 

1604 

.40354 

23 

2191 

-3.5779 

1718 

.19417 

24 

2195 

-3.746 



25 

2197 

-3.89 


Input Screen No, 6 



DATA FOR FREE SURFACE NODES 
Enter: Potential values for mode number 3 


Usa tha ENTER or DOWN ARROW kay to antar aach valua from tha fila. 
U aa tha UP ARROW k.y to novo back up through tha data. 



Noda # 

Fotantial 


Noda # 

Fotantial 

1 

79 

0 

13 

1821 

1.2743 

2 

248 

-.42058 

14 

1900 

1.1995 

3 

416 

-.76613 

15 

1967 

.8525 

4 

582 

-.97616 

16 

2069 

.33672 

5 

745 

-1.0243 

17 

2069 

-.28961 

€ 

904 

-.90966 

18 

2106 

-.97582 

7 

1058 

-.65127 

19 

2135 

-1.6773 

8 

1206 

-.28438 

20 

2157 

-2.3587 

9 

1347 

.14303 

21 

2173 

-2.9891 

10 

1480 

.57584 

22 

2184 

-3.5445 

11 

1604 

.95534 

23 

2191 

-4.0134 

12 

1718 

1.2164 

24 

2195 

-4.4017 




25 

2197 

-4.7374 


V 


Input Screen No. 7 


r 


DATA FOR TANK NALL NODES 





Enter : 

: Node #, and 

R, Z coordinates of 

node, 

starting at R= 

! 


Uaa tha 

ENTER or DOWN 

ARROW 

kay to antar 

aach valua from tha 

a 

H 

%4 


Uaa tha 

UP ARROW kay to mova 

back up through tha 

data. 



Mod* f 

R 

Z 

Mod* # 

R 

z 

1 

1 

0 

-1 

13 

1743 

.9966 

-.0821 

2 

170 

.124 

-.993 

14 

1834 

.9995 

.0322 

3 

338 

.2461 

-.9705 

15 

1912 

.9909 

.1343 

4 

504 

.3644 

-.933 

16 

1978 

.9744 

.2247 

5 

667 

.4772 

-.8811 

17 

2033 

.9526 

.3041 

6 

826 

.5826 

-.8154 

16 

2078 

.9276 

.3735 

7 

980 

.679 

-.7372 

19 

2114 

.9009 

.4341 

8 

1128 

.7649 

-.6475 

20 

2142 

.8735 

.4868 

9 

1269 

.8389 

-.5479 

21 

2163 

.8464 

.5325 

10 

1402 

.9 

-.4398 

22 

2178 

.8201 

.5722 

11 

1526 

.947 

-.3249 

23 

2188 

.7949 

.6068 

12 

1640 

.9794 

-.205 

24 

2194 

.771 

.6367 

V 




25 

2197 

.7487 

.6629 


Input Screen No. 8 
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U) K> 


r 


DATA FOR TANK WALL NODES 
Enter: Potential values for mode number 1 


Oaa tha ENTER or DOWN ARROW kay to antar aach valua from tha fila. 
Osa tha UP ARROW kay to mova back up through t ha data. 



Noda # 

Potantlal 


Noda f 

Potantlal 

1 

1 

0 

13 

1743 

.8799 

2 

170 

.05298 

14 

1834 

.9823 

3 

338 

.1089 

15 

1912 

1.085 

4 

504 

.169 

16 

1978 

1.186 

5 

667 

.233 

17 

2033 

1.282 

6 

826 

.3006 

18 

2078 

1.371 

7 

980 

.3716 

19 

2114 

1.451 

8 

1128 

.4462 

20 

2142 

1.522 

9 

1269 

.5243 

21 

2163 

1.582 

10 

1402 

.6063 

22 

2178 

1.632 

11 

1526 

.6926 

23 

2188 

1.672 

12 

1640 

.784 

24 

2194 

1.702 




25 

2197 

1.727 


v 

Input Screen No. 9 

— \ 

DATA FOR TANK WALL NODES 
Enter: Potential values for mode number 2 


Uaa tha ENTER or DOWN ARROW kay to antar a ach valua from tha fila. 
Uaa tha UP ARROW kay to mova back up through tha data. 



Noda # 

Potantlal 


Noda # 

Potantlal 

1 

1 

0 

13 

1743 

.06863 

2 

170 

.06443 

14 

1834 

-.1741 

3 

338 

.1303 

15 

1912 

-.4842 

4 

504 

.1966 

16 

1978 

-.8512 

5 

667 

.2603 

17 

2033 

-1.259 

6 

826 

.3177 

18 

2078 

-1.688 

7 

980 

. 3646 

19 

2114 

-2.117 

8 

1128 

.3964 

20 

2142 

-2.527 

9 

1269 

.4074 

21 

2163 

-2.902 

10 

1402 

.3904 

22 

2178 

-3.227 

11 

1526 

.3361 

23 

2188 

-3.497 

12 

1640 

.2328 

24 

2194 

-3.712 




25 

2197 

-3.89 




J 


Input Screen No. 10 
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DATA FOR TANK WALL NODES 
Enter: Potential values for mode number 3 

Usa tha ENTER or DOWN ARROW kay to antar aach valua from tha fila. 
Usa tha UP ARROW k«y to iDova back up through tha data. 



Noda # 

Potantial 


Noda # 

Potanti; 

1 

1 

0 

13 

1743 

.738 

2 

170 

-.06434 

14 

1834 

.8465 

3 

338 

-.1249 

15 

1912 

.8167 

4 

504 

-.1757 

16 

1978 

.6092 

5 

667 

-.2083 

17 

2033 

.2139 

6 

826 

-.2142 

16 

2078 

-.3489 

7 

980 

-.1858 

19 

2114 

-1.033 

8 

1126 

-.1169 

20 

2142 

-1.78 

9 

1269 

-.00369 

21 

2163 

-2.53 

10 

1402 

.1533 

22 

2178 

-3.226 

11 

1526 

.3457 

23 

2188 

-3.829 

12 

1640 

.5533 

24 

2194 

-4.323 




25 

2197 

-4.737 


Input Screen No. 11 


r > 

INPUT DATA 

Bond Numbar rill Laval % Cont. Ang Surfaca Maaa Distribution Function 
1.0 50.0 5.0 1.100 - (S/Smax) A 2 . 50 

PARTICIPATION FACTORS OF TRIAL MODES 


B (1) ■ 1.00000 


B (2) m\ 0.01013 


B (3) - 0.00026 


NON-DIMENSIONAL PARAMETERS OF THE PENDULUM MODEL 

(Liq Vol)/{7t* Ro A 3) Sloah M/Liquid M Band. L/Ro 


Ho/Ro 


CG Loc./Ro 


Spring/ (a * Ro A 2) Fraq A 2/ ( (l+Bo)o/ (d * Ro A 3)) 


Do you want tha abova data printad out (Y or N) ? 
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LISTING OF 
COMPUTER PROGRAM 
"GENSURF" 
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ooooooooooooooooooooooooooooo 


c ****************************** PROGRAM ★★*★*★*★********★*★**★ -0001 
c ****************************** G E N S U R F ********************** -0002 


c -0003 

c PURPOSE: Generate the shape of an equilibrium free surface of a -0004 

c liquid at low gravity in sperical and cylindrical -0005 

c containers. -0006 

c -0007 

c USAGE: This program is based heavily on the DYSIM program for -0008 

c system simulation. This routine employs a fourth-order -0009 

c Runge-Kutta algorithm to integrate a system of first-order -0010 

c differential equations. -0011 

c ( ( -0012 

c A file which holds the solution parameters and ODE initial -0013 

c conditions is read. The program requests a value from the -0014 

c user for RLMB (lambda), the pressure jumnp parameter. -0015 

The program solves the surface equations and determines -0016 

the wall contact angle. The user changes the value for -0017 

RLMB until the desired value of wall contact angle is -0018 

achieved. It should be noted that the solution is highly -0019 

sensitive to the value of RLMB. The user should use small -0020 
values (RLMB=0 . 2 ) initially and iterate from there. -0021 

-0022 

A version of this program which employs a Newton-Raphson -0023 

technique for converging on a value of RLMB was prepared. -0024 
This version is highly problem dependent and is not -0025 

available for general use. -0026 

-0027 

MAIN VARIABLES: -0028 

-0029 

t(l) - Length along the surface of the liquid -0030 

t(2) - Upper limit for surface length (prevents runaway solution) -0031 

t ( 3 ) - Integration step size along the surface -0032 

t ( 4 ) - Interval at which to record the solution parameters -0033 

x(l) - Height of the free surface, Z -0034 

x(2) - First derivative of surface height, dZ/ds -0035 

x(3) - Radius of the free surface from tank centerline, R -0036 

x(4) - First derivative of radius, dR/ds -0037 

n - Order of system (always set to n-4) -0038 

m - Number of outputs (set to m«8) -0039 

-0040 

NOTE: The values for surface length, height, and radius are -0041 

normalized with respect to the tank wall radius. -0042 

This routine uses no set of units. -0043 

-0044 

c -0045 

c INPUTS (read from the file called GENSURF . INP ) : -0046 

c ( ( -0047 

c tinit - Initial value of surface length -0048 

c fintim - Final value of surface length -0049 

c step - Integration step size along the surface -0050 

c prtstp - Surface length interval for saving solution parameters -0051 

c issflg - Variable step size flag -0052 

c dtllim - Lower limit for step size -0053 

c dtulim - Upper limit for step size -0054 

c errul - Upper limit on error criteria for variable step size -0055 

c errll - Lower limit on error criteria for variable step size -0056 

c iplot - Flag for line printer plot -0057 

c zO - Initial condition for surface height -0058 

c dzdsO - Initial condition for surface height derivative -0059 

c rO - Initial condition for surface radius -0060 

c drdsO - Initial conditino for surface radius derivative -0061 

c bond - Bond Number -0062 

c thetO - Contact angle at wall -0063 
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c iclspr - Flag for tank type = 0 for spherical tank 
c . ne. 0 for cylindrical tank 

c 

c INPUT from the keyboard: 


c 


c 

c 

rlmb 

— 

Pressure jump parameter 

c 

c 

Suggested 

Values 

for Inputs 

c 

c 

tinit 

_ 

0.0 

Solution should always start at s=0.0 

c 

f intim 

- 

5.0 

Set to a large value, prevents runaway solution 

c 

step 

- 

0.001 

Suggested value 

c 

prtstp 

- 

0.01 

Suggested value 

c 

issflg 

- 

0 

Solution should use constant step (issflg - 0) 

c 




Performance is erratic for variable step (issflg=l) 

c 

dtllim 

- 

.001 

Dummy variable for issflg*0 

c 

dtulim 

- 

.05 

Dummy variable for issflg=0 

c 

errul 

- 

.005 

Dummy variable for issflg=0 

c 

errll 

- 

.001 

Dummy variable for issflg*0 

c 

iplot 

- 

0 

0 for no line printer plot 

c 

zO 


1 

o 

o 

<Tt 

1 for line printer plot at end of output 

c 

— 


c 

dzdsO 

- 

0 

i 

c 

rO 

- 

0 

i 

c 

drdsO 

- 

1 

| These values are for 75% full spherical tank 

c 

bond 

- 

1 

1 

c 

thetO 

- 

5 

1 

c 

c 

iclspr 

rlmb 

— 

0 

2.9806 

1 


c 

c 

c FILES : 
c 

c 1 GENSURF . OUT - Holds ASCII output for line printer 

c 2 GENSURF. ASC - High precision output table use by other routines 

c 3 GENSURF. INP - Parameter input file 

c 

c 

c 

implicit double precision (a-h,o-z) 
c 

common /sim/ x(20) f dx(20) f t(4), y(20) f n, m , iplot, issflg, 

1 dtulim, dtllirn, tprt 

common /param/ bond, rlmb, zO, thetO, pi, iqflg, iclspr 
common/ t save/ tinit 

common /save/ xsavem(20) ,xsave (20) , tsavem,tsave 
common /error/ errmax, ierrmx, ierrf 1, errll, errul, xtrlim 
c 

c 

c Order of system and number of outputs 
c Set the value of pi for degree-radian conversions 



c 

n=4 

m*8 

pi = 4.d0 * atan(l.dO) 

c 

c Open all files here 

c 

open (unit=l, file*' GENSURF . out' , status*' unknown' ) 
open (unit-2 f file*' GENSURF . asc' , status*' unknown' ) 
open (unit=3, file*' GENSURF . inp' , status*' old' , mode*' read' ) 


-0064 

-0065 

-0066 

-0067 

-0068 

-0069 

-0070 

-0071 

-0072 

-0073 

-0074 

-0075 

-0076 

-0077 

-0078 

-0079 

-0080 

-0081 

-0082 

-0083 

-0084 

-0085 

-0086 

-0087 

-0088 

-0089 

-0090 

-0091 

-0092 

-0093 

-0094 

-0095 

-0096 

-0097 

-0098 

-0099 

-0100 

-0101 

-0102 

-0103 

-0104 

-0105 

-0106 

-0107 

-0108 

-0109 

-0110 

-0111 

-0112 

-0113 

-0114 

-0115 

-0116 

-0117 

-0118 

-0119 

-0120 

-0121 

-0122 

-0123 

-0124 

-0125 

-0126 
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c -0127 

c Get all of the input parameters -0128 

c -0129 

read(3,*) tinit -0130 

read(3,*) fintim -0131 

read(3,*) step -0132 

read(3,*) prtstp -0133 

read(3,*) issflg -0134 

read<3,*) dtllim -0135 

read(3,*) dtulim -0136 

read(3,*) errul -0137 

read (3, * ) errll -0138 

read(3,*) iplot -0139 

read (3, *) zO -0140 

read(3,*) dzdsO -0141 

read(3,*) rO -0142 

read(3,*) drdsO -0143 

read(3,*) bond -0144 

read (3, *) thetO -0145 

read(3,*) iclspr -0146 

c -0147 

c Get the jump parameter from the user -0148 

c -0149 

print*,' Enter the guessed value of pressure jump condition' -0150 

print*,' for the desired contact angle=' , thetO -0151 

read(5,*) rlmb -0152 

c -0153 

c Transfer inputs to program variables and other setup stuff -0154 

c -0155 

t(l) -tinit -0156 

t (2 ) = fintim -0157 

t (3) = step -0158 

t ( 4 ) = prtstp -0159 

x { 1 ) - zO -0160 

x(2) = dzdsO -0161 

x (3) - rO -0162 

x ( 4) * drdsO -0163 

thetO = thetO * pi/180. -0164 

xtrlim - (t (2)-t (1) )/t (3) * 10. -0165 

c -0166 

c Record input variables for the user -0167 

c -0168 

write (1,16010) -0169 

write (1, 16020) tinit, fintim, step, prtstp -0170 

write (1, 16030) issflg, errul, errll, dtulim, dtllim, xtrlim -0171 

write (1, 16040) iplot -0172 

write (1, 16050) x(l), x(2), x(3), x(4) -0173 

write (1, 16060) bond, thet0*180 . /pi, rlmb -0174 

if (iclspr .eq. 0) then -0175 

write (1, 16070) iclspr -0176 

else -0177 

write (1, 16080) iclspr -0178 

endif -0179 

-0180 

-0181 

c -0182 

c Begin the solution -0183 

c -0184 

call dysim -0185 

c -0186 

c Report the results to the user -0187 

c -0188 

if (iclspr .eq. 0) then -0189 
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angl = 

-pi/2. 




-0190 

ang2 * 

pi/2. 




-0191 

if (dx ( 1 ) .gt. 0.) angl - atan (-dx (3) /dx ( 1) ) 



-0192 

if (x(l) .It. 1.) ang2 = asin(x(l)) 




-0193 

angle 

= ang2 - angl 




-0194 

else 





-0195 

angle 

= pi/2 . 




-0196 

if (abs(dx(l)) .le. 1.) angle = acos(dx(l)) 



-0197 

endif 





-0198 

c 





-0199 

print*, ' 

tried for contact angle = ',thet0* 

180 . /pi 



-0200 

print*, ' 

found contact angle = ', angle * 180. /pi 



-0201 

print*, ' 

rlmb = ' , rlmb 




-0202 

write (1, 16999) rlmb, angle*180 . /pi 




-0203 

c 





-0204 

stop 





-0205 

c 





-0206 

16010 format (lhl/' Liquid Free Furface' ) 



-0207 

16020 format (///5x, ' Starting Value of Surf. Length [t(l)] 

= 

t 

r 

-0208 

1 

t48, lpel2 . 5 , 9 sec' , / 




-0209 

1 

5x, 'Final Value of Surf. Length 

[t (2) ] 

= 

t 

r 

-0210 

1 

t48, lpel2 . 5, ' sec' , / 




-0211 

1 

5x, 'Step Size 

[t (3) ] 


r 

r 

-0212 

1 

t48 , lpel2 . 5, ' sec' , / 




-0213 

2 

5x, ' Print/Plot interval 

[t (4) ] 

= 

r 

9 

-0214 

3 

t48, lpel2 . 5, ' sec' , /) 




-0215 

16030 format (/ 

5x, 'Step Size Control Flag 

(issflg) 


' , t48, i2, / 

-0216 

1 

5x, ' Upper Limit on Error Estimate 

(errul) 

= 

9 

9 

-0217 

1 

t 48, lpel2 . 5, / 




-0218 

2 

5x, ' Lower Limit on Error Estimate 

(errll) 

= 

9 

r 

-0219 

1 

t48 , lpel2 . 5, / 




-0220 

1 

5x, 'Upper Limit on Step Size 

(dtulim) 

= 

9 

9 

-0221 

1 

t 48 , lpel2 . 5, ' sec' , / 




-0222 

1 

5x, 'Lower Limit on Step Size 

(dtllim) 

- 

9 

9 

-0223 

1 

t48 , lpel2 . 5, ' sec',/ 




-0224 

1 

5x, 'Maximum Number of Steps 

(xtrlim) 

* 

9 

9 

-0225 

1 

t48 , lpel2 . 5, / ) 




-0226 

16040 format (/ 

5x, ' Printer/Plot flag 

(iplot) 

* 

' ,t48,i2 f /) 

-0227 

16050 format (/ 

5x, ' Initial Height 

(x(l)] 

= 

9 

9 

-0228 

1 

t48, lpel2 . 5, ' ',/ 




-0229 

1 

5x, ' Initial dZds 

[x (2) ] 

= 

9 

9 

-0230 

1 

t48, lpel2 . 5, ' ',/ 




-0231 

1 

5x, 'Initial Radius 

[x (3) ] 

«= 

9 

9 

-0232 

1 

t 48, lpel2 . 5, ' ',/ 




-0233 

1 

5x, ' Initial dRds 

[x(4) ] 

= 

9 

9 

-0234 

1 

t48 , lpe!2 . 5, ' ',/) 




-0235 

16060 format (/ 

5x,'Bond Number 

(bond) 

* 

9 

9 

-0236 

a 

t48 , lpel2 .5, ' ',/ 




-0237 

1 

5x, 'Desired Contact Angle 

(thet) 


9 

9 

-0238 

a 

t48, lpel2 .5, / 




-0239 

1 

5x, 'Pressure Jump condition 

(rlmb) 


9 

9 

-0240 

a 

t48, lpel2 . 5, / ) 




-0241 

16070 format (/ 

5x, 'Tank shape flag 

(iclspr) 

* 

9 

9 

-0242 

1 

t48,i3,' Spherical Tank' 

) 



-0243 

16080 format {/ 

5x, 'Tank shape flag 

(iclspr) 

= 

9 

9 

-0244 

1 

t48,i3,' Cylindrical Tank') 



-0245 

16999 format (///' For Pressure Jump Parameter (rlmb) = ' , 

lpe!2 .5/ 

-0246 

1 

' Found Contact Angle = ',lpel2.5) 



-0247 

end 





-0248 
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c *********************** SUBROUT INE ************************ —0249 
0 ★*★*★★★★*★★★★★★★★*★★*★**★* * D E R F U N ★★*★*★★★★★*★*★***★**★*★★ —0250 

subroutine derfun -0251 

c “0252 

c Defines the derivative functions to be integrated -0253 

c -0254 

c These functions are the first order derivatives that define the -0255 

c system of model equations to be solved. -0256 

c -0257 

c -0258 

implicit double precision (a-h,o-z) -0259 

c -0260 

common /sim/ x(20), dx(20), t(4), y(20), n, m , iplot, issflg, -0261 

1 dtulim, dtllim, tprt -0262 

common /param/ bond, rlmb, zO, thetO, pi, iqflg, iclspr -0263 

common/tsave/ tinit -0264 

common /save/ xsavem(20) , xsave (20) , tsavem, tsave -0265 

common /error/ errmax, ierrmx, ierrf 1, err 11, errul, xtrlim -0266 

c -0267 

c -0268 

c -0269 

c Derivative functions -0270 

c -0271 

dx (3) = x ( 4 ) -0272 

dx ( 1 ) = x (2 ) -0273 

if (t(l) .eq. 0.) then -0274 

dx (2 ) = 0.5 * rlmb -0275 

else -0276 

dx (2 ) = ( bond* (x(l) -zO) + rlmb)*dx(3) - dx (1) *dx (3) /x (3) -0277 

endif -0278 

if (x (3 ) .ne. 0.) then -0279 

dx ( 4) = -(bond*(x(l)-z0)+rlmb)*dx(l) + dx (1) *dx (1) /x (3) -0280 

else -0281 

dx (4) = 0. -0282 

endif -0283 

c -0284 

c -0285 

c Check radius to determine if we are at tank wall -0286 

c -0287 

iqflg = 0 -0288 

wllchk = sqrt (x (1) *x ( 1) + x(3)*x(3)) -0289 

if (iclspr .ne. 0) wllchk = x(3) -0290 

if (wllchk .ge. 1. .or. x(3) . le . 0.) iqflg = 1 -0291 

if (iclspr .ne. 0 .and. dx(3) . le . 0.) iqflg = 1 -0292 

c -0293 

c Enter outputs to be printed = y(l) through y(m) -0294 

c -0295 

c Originally: -0296 

c -0297 

c y(l) = Surface height, Z -0298 

c y(2) = Surface height first derivative, dZ/ds -0299 

c y(3) = surface height second derivative, d2Z/ds2 -0300 

c y (4) = Surface radius, R -0301 

c y(5) = Surface radius first derivative, dR/ds -0302 

c y(6) ~ surface radius second derivative, d2R/ds2 -0303 

c y(7) = local angle of surface with the 'horizontal' -0304 

c y (8) = wall contact angle (valid only at last point in solution) -0305 

c -0306 

y(l) = x (1) -0307 

y (2) = x (2 ) -0308 

y (3) - dx (2) -0309 

y (4) = x (3) -0310 

y (5) - x ( 4 ) -0311 
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y ( 6) = dx (4) -0312 

y (7) = 90. -0313 

if (dx(3) .ne. 0.) then -0314 

y (7 ) = at an (dx (1) /dx (3) ) * 180. /pi -0315 

endif -0316 

if (iclspr .eq. 0) then -0317 

angl = -pi/2. -0318 

ang2 = pi/2. -0319 

if (dx(l) .gt. 0.) angl = atan (-dx (3) /dx (1) ) -0320 

if ( x ( 1 ) .It. 1.) ang2 = asin(x(l)) -0321 

y ( 8 ) - (ang2 - angl) * 180. /pi -0322 

else -0323 

y ( 8) = pi/2 . -0324 

if <abs(dx(l)) .le. 1.) y(8) * acos(dx(l)) -0325 

c if (dx ( 1) -ne. 0.) y(8) - atan (dx (3) /dx (1) ) -0326 

y (8) = y (8 ) * 180. /pi -0327 

endif -0328 

c -0329 

return -0330 

end -0331 
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c ** **★*★★**★*★****★*★** * SUBROUT INE ************************ 
c *************************** D Y S I M ************************ 

c 

c Controls the solution process 
c 

c > This routine should not be changed by the user < ! ! ! ! ! 

c 

subroutine dysim 


c 

c 


c 


c 


c 

c 

c 

c 


implicit double precision (a-h,o-z) 


common /sim/ x{20), dx(20), t(4) f y(20), n, m , iplot, issflg, 
1 dtulim, dtllim, tprt 

common /param/ bond, rlmb, zO, thetO, pi, iqflg, iclspr 
common/tsave/ tinit 

common /save/ xsavem(20) , xsave (20) , tsavem, tsave 
common /error/ errmax, ierrmx, ierrf 1, errll, errul, xtrlim 


dimension tl(1500), yp(20,1500), ymin(20), ymax(20) 
integer*2 iy (20 ) , plot ( 101) , blank, star (20) 


data iy, blank /20*'y(',' '/ 
data star / r 1 9 ,'2', '3', '4', '5', '6', '7', '8' 
1 'B' , ' C' , ' D' , 'E' , 'F' , 'G' , 'H' , ' I' 


9' ,'A', 
J'/K'/ 


Setup for solution — Print initial conditions 


call derfun 


tprt = t (1) 

if ( (t(2)-t(l))/t(4) .le. 1500.) go to 10 
write (1,16010) (iy(i),i,i=l,m) 

iplot=0 
10 lp-1 

tl(lp) - t(l) 
do 20 i=l,m 

yp(irip) = y(i) 

20 continue 

write (1,16001) (iy(i),i,i=l,m) 

nflin = (m-10)/10 
nllin « mod (m, 10) 
ml = 10 

if (m .le. 10) then 
ml - m 

write (1, 16020) t(l), (y (i) , i=l,ml) 
else 

write (1,16020) t(l), (y (i) , i=l,ml) 
if (nflin .eq. 0) then 

write (1,16030) (y(i), i=ll,m) 

else 

if (nllin .eq. 0) then 

write (1,16040) (y(i),i=ll,m) 

else 

write (1,16050) (y(i),i=ll,m) 

endif 
endif 
endif 

do 30 i=l,m 

ymin(i) = 0.9999*y(i) 
ymax(i) = 1.0001*y(i) 

30 continue 
k = 0 


-0332 

-0333 

-0334 

-0335 

-0336 

-0337 

-0338 

-0339 

-0340 

-0341 

-0342 

-0343 

-0344 

-0345 

-0346 

-0347 

-0348 

-0349 

-0350 

-0351 

-0352 

-0353 

-0354 

-0355 

-0356 

-0357 

-0358 

-0359 

-0360 

-0361 

-0362 

-0363 

-0364 

-0365 

-0366 

-0367 

-0368 

-0369 

-0370 

-0371 

-0372 

-0373 

-0374 

-0375 

-0376 

-0377 

-0378 

-0379 

-0380 

-0381 

-0382 

-0383 

-0384 

-0385 

-0386 

-0387 

-0388 

-0389 

-0390 

-0391 

-0392 

-0393 

-0394 
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C -0395 

c Write initial conditions to high precision output file -0396 

c -0397 

write (2 , 16060) t(l), (y(i),i=l,7) -0398 

c -0399 

c Solution loop -0400 

c -0401 

100 continue -0402 

k * k + 1 -0403 

if (k .gt. xtrlim) go to 999 -0404 

ierrfl = 0 -0405 

call rkint -0406 

if (ierrfl .ne . 0) go to 100 -0407 

do 110 i=l , m -0408 

if (y ( i ) .It. ymin(i)) ymin(i) * y(i) -0409 

if (y ( i ) .gt. ymax(i)) ymax(i) = y(i) -0410 

110 continue -0411 

if ( t(l) .ge. tprt .or. iqflg .eq. 1) then -0412 

tprt = tprt + t ( 4) -0413 

nflin = (m-10)/10 -0414 

nllin = mod (m, 10) -0415 

ml = 10 -0416 

if (m .le. 10) then -0417 

ml - m -0418 

write (1,16020) t (1) , (y (i) , i=l,ml) -0 419 

else -0420 

write (1, 16020) t(l), (y (i) , i-1, ml) -0421 

if (nflin .eq. 0) then -0422 

write ( 1, 16030) (y (i) , i=ll,m) -0423 

else -0424 

if (nllin .eq. 0) then -0425 

write ( 1, 16040) (y (i) , i^ll, m) -0426 

else -0427 

write {If 16050) (y (i) , i=ll, m) -0428 

endif -0429 

endif -0430 

endif -0431 

c -0432 

c Write to unformatted file -0433 

c -0434 

write (2, 16060) t (1) , <y(i) ,i=l, 7) -0435 

c -0436 

c Save stuff for plotting if needed -0437 

c -0438 

if (iplot .eq. 1) then -0439 

lp *= Ip + 1 -0440 

tl(lp) = t(l) -0441 

do 120 i«l,m -0442 

yp (i, Ip) * y (i) -0443 

120 continue -0444 

endif -0445 

endif -0446 

c -0447 

c Test for end -0448 

c -0449 

if (iqflg .ne. 1) then -0450 

if (t { 1 ) .It. t (2 ) ) then -0451 

if ( (t (2) -t ( 1) ) .gt. t (3) ) go to 100 -0452 

t (3) = t (2) - t(l) -0453 

tprt = t ( 4) -0454 

go to 100 -0455 

endif -0456 

endif -0457 
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c -0458 

c Produce the line printer plot -0459 

c -0460 

if (iplot .eq. 0) go to 300 -0461 

write (1, 16070) (i, ymin (i) , i f ymax (i) , i=l, m) -0462 

write (1,16080) -0463 

do 230 1=1, lp -0464 

do 210 j=l, 101 -0465 

210 plot(j) = blank -0466 

do 220 i=l,m -0467 

jp = 1 -0468 

denom = ymax(i) - ymin(i) -0469 

if (denom .ne. 0.) then -0470 

jP = (yp (i, 1) -ymin (i) ) /denom * 100 + 1 -0471 

endif -0472 

plot ( jp) = star (i) -0473 

220 continue -0474 

write (1, 16090) tl(l),plot -0475 

230 continue -0476 

c -0477 

300 continue -0478 

c -0479 

return -0480 

c -0481 

c -0482 

c Error conditions -0483 

c -0484 

999 continue -0485 

write (1,169 99) float (k) , xtrlim -04 8 6 

print 16999, float (k) , xtrlim -0487 

stop -0488 

c -0489 

16001 f ormat ( ' ITank Free Surface Solution'// -0490 

1 5x, 'T (1) ' , 5x, 10 <a2, i2, ' ) ' , 6x) /15x, 10 (a2 , i2, ' ) ' , 6x) ) -0491 

16010 f ormat ( lx,'***** Error Message from DYSIM ****** ,/, -0492 

*' The number of data points to be plotted exceeds the dimensioned -0493 
*number of 1500.'/' Decrease final time or increase print interval. -0494 
*Plot will be supressed.') -0495 

16020 format (lx, lpdlO .3, lplOdll .3) -0496 

16030 format (15x, lplOdll . 3 ) -0497 

16040 format (10 (15x, lplOdll. 3/) , 15x, lplOdll. 3 ) -0498 

16050 format <10 (15x, lplOdll. 3/), 15x, lplOdll. 3 ) -0499 

16060 format (Ip8d20 . 10) -0500 

16070 f ormat (' lPlot of y(i) versus t(l)'// -0501 

* ( 12x, ' YMIN ( ' , i2 , ' ) * ',lpdl0.3,5x,'YMAX(',i2,') =',lpdl0.3)) -0502 

16080 format {/5x, 'time' , tl3, ' 0 . ' ,34x, ' (y (i) -ymin(i) ) / (ymax(i) -ymin (i) ) ' , -0503 

* tll4,'l.'/13x,10(' ! '),'!') -0504 

16090 format (lx, Ipdll . 4, lx, 101a 1) -0505 

16999 format (Ih ********* Program Terminated ********'/ -0506 

1 ' Iteration count = ',lpdl2.4/ -0507 

2 ' Greater than limit of',lpdl2.6) -0508 

c -0509 

end -0510 


B-23 



c *********************** SUBROUT INE ★***★*★★*********★★**★★* -0511 

£★**★***★*★*★***★*★*****★*★* R K I N T ************************ -0512 

c -0513 

c Performs the integration of the derivatives at each step in the -0514 

c solution -0515 

c -0516 

c > This routine should not be changed by the user < ! ! ! ! ! “0517 

c -0518 

c -0519 

subroutine rkint -0520 

c -0521 

implicit double precision (a-h,o-z) -0522 

c -0523 

common /sim/ x(20), dx(20), t(4), y<20), n, m , iplot, issflg, -0524 

1 dtulim, dtllim, tprt ( -0525 

common /param/ bond, rlmb, zO, thetO, pi, iqflg, iclspr -0526 

common /save/ xsavem (20 ) , xsave (20 ) , tsavem, tsave ^ -0527 

common /error/ errmax, ierrmx, ierrf 1, err 11, errul, xtrlim -0528 

c -0529 

dimension a(4,20),c(4) -0530 

c -0531 

c(l) - 0.5 -0532 

c (2 ) - 0.5 -0533 

c (3) = 1.0 -0534 

c (4) = 0.0 -0535 

do 10 i=l, n -0536 

10 xsave (i) * x(i) -0537 

tsave = t(l) -0538 

do 40 k*l, 4 -0539 

call derfun -0540 

do 20 i=l , n -0541 

20 a (k, i) = dx(i) -0542 

do 30 i~l,n -0543 

30 x(i) = xsave(i) + c (k) *a (k, i) *t (3) -0544 

t(l) = tsave + c(k)*t(3) -0545 

40 continue -0546 

errmax = 0. -0547 

do 50 i=l, n -0548 

x(i) = xsave(i) + -0549 

1 t (3) /6.0*(a(l,i)+2.0*a(2,i)+2.0*a(3,i)+a(4,i) ) -0550 

denom = a (2 , i) -a (1, i) -0551 

xnum = a(3,i) - a(2,i) -0552 

iiii = 0 -0553 

if (denom .ne. 0.) then -0554 

err =2.* abs ( (a (3, i) -a (2, i) ) /denom ) -0555 

iiii = i -0556 

elseif (xnum .eq. 0) then -0557 

err = errll*0.998 -0558 

else -0559 

err * errll*1.002 -0560 

endif -0561 

errmax = max (err, errmax) -0562 

if (errmax .eq. err) ierrmx = iiii -0563 

50 continue -0564 

t(l) *= tsave + t(3) -0565 

c -0566 

c check step size -0567 

c -0568 

ierrfl - 0 -0569 

if (issflg .ne. 1) go to 506 -0570 

if (errmax ,gt. errul) then -0571 

t (3) - t (3) * 0.5 -0572 

if (t (3) .le. dtllim) then -0573 
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t(3) = dtllim 


-0574 



endif 


-0575 



t (1) = tsave 


-0576 



do 505 i=l,n 


-0577 


505 

x ( i ) = xsave(i) 


-0578 



ierrfl = 1 


-0579 



go to 510 


-0580 



endif 


-0581 



if (errmax .It. errll) t(3) 

- t<3) * 1.05 

-0582 



if (t(3) .gt. dtulim) t(3) 

= dtulim 

-0583 


506 

continue 


-0584 

c 




-0585 



call derfun 


-0586 

c 




-0587 


510 

continue 


-0588 

c 




-0589 



return 


-0590 



end 


-0591 


B-25 



B-26 



LISTING OF 
COMPUTER CODE 
"LOW-G.BAS" 
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LOW-G . BAS 


This program uses the trial functions from ADINA to compute the low-g slosh natural 
frequency, the slosh mode shape, the slosh force and torque, and the parameters 
of an equivalent mechanical pendulum model. No more than THREE trial modes are 
allowed. The number of nodal points along the surface line THETA = 0 must be 
>= 13 and <- 25. 


Subroutines defined used in this program. 


DECLARE SUB AIntegrals (F!(), APartial!(), Coord! (), WH!(), Result! ()) 
DECLARE SUB APart (Coord! (), Kappa2 ! ( ) , WD2 ! <) , WH!(), ResultlO) 
DECLARE SUB BIntegrals (Coord! (), Func!(), WH!(), ResultlO) 

DECLARE SUB CompDerivs (Coord! {) , Value! () , Deriv! {) ) 

DECLARE SUB Deriv (Code%, Coord! (), FuncDat ! (), FuncDerivI () ) 

DECLARE SUB InputScreen (FileNo%, Code% f Temp$(), TempNoO) 

DECLARE SUB LiqVol (Surf!<), Wall!(), Ansi!, Ans2 ! , DZ ! , DeltaElev!, 
SurfElev ! ) 


DECLARE 
DECLARE 
DECLARE 
DECLARE 
DECLARE 
DECLARE 
DECLARE 
DECLARE 
DECLARE 
Ans2 ! ) 


FUNCTION PadZero$ (A!, NoDecs ! , NoBef oreDec ! ) 

SUB ParamModel (MOF ! , F2 ! f SM! , SM2 ! , PL!, SC!, HO!) 

SUB PotNorm (SurfP!(), WallP ! () ) 

SUB Printout (F2 ! , F!(), SM2 ! , PL!, HO!, CG!, SC!, BI!(), PC%) 

SUB PrintOut2 (S!<), WH ! ( ) , BI!(), F2 ! , F ! ( ) ) 

SUB OutputScreen (F2 ! , SM2 ! , PL!, HO!, CG!, SC!, BI!(), PCode%) 

SUB Rad (Coord! (), Temp!()) 

SUB Root (Guess!, DT!, Cl!(), C2!(), Answer!) 

SUB SloshM (S ! ( ) , W!(), PS!(), PW!(), TF!(), BI!(), SlsFrq!, Ansi!, 


DECLARE SUB Text Ini <T$, Max%, NumOnly%, CapsOn%, ExitCode%, Colr%) 
DECLARE SUB Vector (F2 ! , Cl!<), C2!(), Result! ()) 

DECLARE SUB WaveHite (Coord! (), NodeDatS!(), TempNo!()) 


The following routines are assembly-language routines from the QUICKPAK library. 
Both are used only in the screen input routine f, TextInl" . "Peekl%" reads a byte in 
memory; its use here is to determine what kind of monitor is being used, so it can 
be deleted without problem. "QPrint" just prints a text string on the monitor 
screen rapidly, so it can be replaced if needed by a combination of regular LOCATE 
statements and PRINT statements. 


DECLARE FUNCTION Peekl% (segment. Address) 
DECLARE SUB QPrint (X$, Colr%, Page%) 


DIM ModeS (4, 25), ModeW(4, 25), WavHite(3, 25) 

DIM CoordS (7, 25), CoordW(3, 25), RadCurv2(25) 

DIM TempStor$ (10, 25) 

DIM Freq (4) 

DIM WavDerivl (3, 25), WavDeriv2 (3, 25) 

DIM Al(3, 25), A ( 3, 3), B<4, 4), Vec(3) 

COMMON NoModes%, ModeNo%, NoNodes%, ConstMass, ExpMass, BondNo 
COMMON Ang, FillLevel, LiqVolume, LastDZsurf, DeltaElev, SurfElev 

CLS ' clear the screen for the data input routines 

' ================ SEGMENT 1 SEGMENT 1 SEGMENT 1 ================ 

' The code in Segment 1 does the following: (1) puts the name of the program on the 

' screen, (2) asks for the number of trial modes and the number of nodal points 
' for each mode along the free surface and the wall, and (3) whether the data will 
' be entered from the keyboard or read in from disk files named COORD .DAT and 


PRECEDING PAGE BLANX NOT FILMED 




' MODE.DAT. The code uses some subroutines from QUICKPAK to draw boxes around the 
' instruction and input text; the subroutine can be deleted if the QUICKPAK 
' library is not available. 

Start : 

' Print the program title in a box 

CALL Box < 1, 5, 6, 75, 2, 12, -1) 

COLOR 14: LOCATE 3, 14 

PRINT "MINIMIZE INTEGRAL TO FIND THE LOW-G SLOSH POTENTIAL" 

LOCATE 4, 17 

PRINT "BY USING THE STRUCTURAL CODE TRIAL FUNCTIONS" 

' Draw boxes and give instructions for entering the input data: 

B$ - SPACES (75) use BS to clear parts of the screen later 

COLOR 7: LOCATE 7, 25: PRINT "Hit RETURN after each entry" 

Responses - "N" 

WHILE Responses = "N" ' repeat instructions until input data is okay 


CALL Box (8, 

17, 

10, 22, 1, 

14, -1) 

LOCATE 9 , 6 , 

1, 

5, 7: PRINT 

"Bond No. ■= " 

CALL Box (8, 

36, 

10, 39, 1, 

14, -1) 


LOCATE 9, 24: PRINT "Filling % = " 

CALL Box ( 8 , 73, 10, 76, 1, 14, -1) 

LOCATE 9, 43: PRINT "Contact angle at wall (deg) = " 

CALL Box (11, 32, 13, 37, 1, 14, -1) 

LOCATE 12, 6: PRINT "Mass function: Exponent = " 

CALL Box(ll, 54, 13, 60, 1, 14, -1) 

LOCATE 12, 42: PRINT "Constant - " 

' Now get the input: 

LOCATE 9 , 18: INPUT "", BondNo ' Bond number 

LOCATE 9, 37: INPUT FillLevel ' tank fill level 0 - 100% 

LOCATE 9 , 74: INPUT Ang ' contact angle at the wall (deg) 

LOCATE 12, 33: INPUT ExpMass ' constant in the mass distribution 

LOCATE 12, 55: INPUT "" , ConstMass ' exponent in the mass distribution 

' Is the data okay? 

LOCATE 14, 10 

INPUT "Are all of the input numbers okay (Y,N) ? ", Response? 
Response? * UCASE? (Response?) ' convert to uppercase 

IF Response? * "N" THEN 

LOCATE 8, 1 ' clear out the data 

FOR I = 1 TO 6 
PRINT B? 

NEXT I 
END IF 

WEND ' end of the "WHILE-WEND" input loop 

LOCATE 14, 1: PRINT B? ' clear out the line of text 

' Get the number "NoModes%" of trial modes and the number M NoNodes%" of 
' nodal values per mode along THETA=0 from the structural simulation, and 
' the trial natural frequencies: 

Response? * "N" 

WHILE Response? = "N" ' do until the data is okay 

CALL Box{14, 52, 16, 54, 1, 14, -1) 

LOCATE 15, 6 

PRINT "How many trial modes do you want to consider? " 

CALL Box (16, 64, 18, 67, 1, 14, -1) 

LOCATE 17, 6 

PRINT "How many nodal values are there at the surface and wall? " 

Mode Again: 

LOCATE 15, S3: INPUT "", NoModes% 

IF NoModes% > 3 THEN CALL Chime (6) : GOTO ModeAgain 
Node Again: 

LOCATE 17, 65: INPUT "", NoNodes% 



IF NoNodes% < 13 OR NoNodes% > 25 THEN CALL Chime (6) : GOTO NodeAgain 
LOCATE 20, 6: PRINT "Eigenvalue of trial modes — " 

FOR NF% * 1 TO NoModes% 

IColumn% = 23 + 11 * NF% 

IColumn2% = IColuinn% + 10 

CALL Box (19, IColumnl, 21, IColumn2%, 1, 14, -1) 

LOCATE 20, IColumn% + 1: INPUT Freq(NF%) 

NEXT NF% 

LOCATE 22, 10 

INPUT "Are these last input numbers okay (Y,N) ? ", Response$ 

Response$ = UCASE$ (Response$) 9 convert to uppercase 

IF Response$ = "N" THEN 
LOCATE 14, 1 

FOR 1 = 1 TO 10 ' clear out the data 

PRINT B$ 

NEXT I 
END IF 

WEND ' end of "WHILE-WEND" loop for modal data 

LOCATE 22, 1: PRINT B$ 

' Find out whether whether the nodal coordinates, trial potential values, 

' etc., will be given from the terminal or from data files: 

FileAgain : 

CALL Box ( 22 , 76, 24, 78, 1, 14, -1) 

LOCATE 23, 6 

PRINT "DO YOU WANT TO INPUT NODAL DATA FROM THE (1) TERMINAL OR (2) FILES ? " 
LOCATE 23, 77: INPUT "", Response$ 

Answer = VAL (Respon se$) ' If Response$=l the data will be typed in 

IF Answer = 1 OR Answer = 2 THEN 
GOTO Cont 
ELSE 

CALL Chime (6) ' input a number <> 1 or 2 or other bad answers 

GOTO FileAgain 
END IF 


' ========= end OF SEGMENT 1 END OF SEGMENT 1 ================== 

' ========== SEGMENT 2 SEGMENT 2 SEGMENT 2 ===================== 

' This segment of code gets (1) coordinates of the equilibrium free surface, 

' (2) the required derivative values of the equilibrium free surface, (3) 

' the coordinates of the finite element nodes along THETA=0, and (4) the 
' values of the trial potentials at these nodes. 

Cont : 

SELECT CASE Answer 

CASE 1 ' Type INPUT at the terminal 

' ================ Option to type in the data input ================== ' 

The following segment gets the typed in data using subroutines. The 

' subroutines use the QUICKPAK routine TEXTIN ! . BAS to allow data 
' correction anywhere on the screen, without including a lot of code. 

' The subroutine can be deleted if QUICKPAK is not available but the 
9 input requests will have to be redone to allow a mistake made in 
' typing in the input to be corrected. The input is saved in files 
' COORD . DAT {for the R, S coordinates and R',Z',R H ,Z" of each nodal 
' point) and M0DE.DAT {for the potential values at each node) . 

' Get data for surface nodes and store in array "CoordS ()." Row 1 
' of CoordS stores the node number, row 2 stores the R coord, row 3 
' stores the S coord. The derivative R' is stored in row 4, Z' is 
' stored in row 5, R n is stored in row 6, and Z" is stored in row 7. 

' The first screen gets the node number and R and Z. The second screen 
' gets R' and Z' . The third screen gets R" and Z". 



CALL Input Screen (O f 1, TempStor$(), CoordS ()) ' get R and S 

CALL InputScreen (0, 2, TempStor$ O , CoordSO) ' R' and Z' 

CALL InputScreen (0, 3, TempStor$<), CoordS ()) ' R" and Z" 

' Get the potential values at each surface node and store in array 
' ModeSO. Row 1 is the node # from CoordSO , row 2 the potential 
' values for the first mode, row 3 the values for the second mode, 

' and so on — Row 1 + 1 is the Mode I data. 

FOR K - 1 TO NoModes% 

ModeNo% = K 

CALL InputScreen (0, 4, TempStor$(), ModeSO) 

NEXT K 

' Get the data for the R, Z coordinates of the wall nodes and store in 
' array CoordWO. Row 1 is the node # data. Row 2 is R, and Row 3 is Z. 

' First, though, empty out the unneeded surface data in TempStor$. 

FOR I * 1 TO 10 
FOR J = 1 TO 25 

TempStor$ (I, J) - "" 

NEXT J 
NEXT I 

CALL Input Screen { 0 , 5, TempStor$(), CoordWO) ' R and Z 

' Get potential values for each wall node point and store in ModeW(). 

' Row 1 is the node #, Row 2 is the first mode, and so on. 

FOR K = 1 TO NoModes% 

ModeNo% = K 

CALL InputScreen (0, 6, TempStor$ {) , ModeW () ) 

NEXT K 

' =========== End of option to type in the input data =========== 

' =========== Start of option to enter data from disk files ========= 

CASE 2 

' Read the surface and wall coordinate data from COORD .DAT file into 
' the array TempStor$ O . Rows 1 - 7 of TempStor$<) will contain the 
' surface node #'s and coord, data, and rows 8-10 will contain the 
' wall node #'s and coord, data. 

OPEN "COORD. DAT” FOR INPUT AS #1 
FOR I = 1 TO 10 

FOR J = 1 TO NoNodes% 

INPUT #1, TempStor$ (I, J) 

NEXT J 
NEXT I 
CLOSE #1 

' Send the surface coordinate data forward for display on the screen 
9 and ask the user if the data is okay. Note that COORD.DAT is 
9 erased and then rewritten after verification. 

CALL InputScreen (1, 1, TempStor$(), CoordSO) 

CALL InputScreen (1, 2, TempStor$ ( ) , CoordSO) 

CALL InputScreen (1, 3, TempStor$ () , CoordSO) 

9 Read the surface and wall modal values form MODE.DAT. The surface 
' values are stored in TempStor${), without deleting the surface 
' node numbers. The wall values are stored temporarily in the array 
' ModeW () and will later be read into TempStor$(). 

OPEN "MODE.DAT” FOR INPUT AS #1 
FOR J = 1 TO 2 * NoModes% 

FOR I * 1 TO NoNodes% 

IF J * 1 THEN TempStor$ * STR$ (CoordS ( 1, I)) ' node # 

IF J <= NoModes% THEN INPUT #1, TempStor$ ( J + 1, I) 

IF J > NoModes% THEN 



JJ = J - NoModes% 

INPUT #1, ModeW ( JJ, I) 

END IF 
NEXT I 
NEXT J 
CLOSE #1 

' Send the surface mode data forward for display on the screen and 
' ask the user if the data is okay. MODE . DAT is erased and then 
' rewritten with the verified data. 

FOR I = 1 TO NoModes% 

ModeNo% = I 

CALL InputScreen (1, 4, TempStor$ ( ) , ModeS ()) 

NEXT I 

9 Now move the wall coordinate data from the last three rows of 
' TempStor${) to the first three rows, so that the data will be 
' displayed correctly on the screen. Note that the last three rows 
' were not overwritten by the surface data input routines above. 

FOR I = 1 TO NoNodes% 

TempStor$ ( 1, I) = TempStor$(8, I) 

TempStor$(2, I) = TempStor$(9 f I) 

TempStor$ (3 , I) = TempStor$ ( 10, I) 

NEXT I 

' Send the wall coordinate data forward for display on the screen and 
' ask the user if the data is okay. COORD . DAT is not erased by this 
' operation, but the verified data overwrites the previous wall data. 
CALL InputScreen {1, 5, TempStor$(), CoordWO) 

' Move the wall coord, data from the array CoordWO to the first row 
9 of TempStor$ O and the wall modal data from the array ModeW () to 
f rows 2, 3, and 4 . 

FOR I - 1 TO NoModes% 

FOR J = 1 TO NoNodes% 

IF I - 1 THEN TempStor$ ( 1, J) - STR$ (CoordW (1, J) ) 
TempStor$ (I + 1, J) = STR$ (ModeW (I, J) ) 

NEXT J 
NEXT I 

9 Send the wall modal data forward for display on the screen and ask 
' the user to verify the data. MODE . DAT is not erased by this 
' operation, but the verified data overwrites the previous wall data. 
FOR K = 1 TO NoModes% 

ModeNo% - K 

CALL InputScreen (1, 6, TempStor$(), ModeWO) 

NEXT K 

' All the data has been displayed, checked, rewritten to COORD . DAT and 
' MODE . DAT, and stored in Coord {) and Mode() arrays. 


' ========= End of option to enter data from disk files === 

CASE ELSE r Just a safety check 

CLS : CALL Chime (6): GOTO Start 

END SELECT ' end of loops to get input data 

' =========== END OF SEGMENT 2 END OF SEGMENT 2 ============ 


/ =========== SEGMENT 3 SEGMENT 3 SEGMENT 3 ============= 

9 This segment does the numerical integrations. First, the trial potentials 
' are normalized to have a maximum magnitude of one. 

CALL PotNorm (ModeS () , ModeWO) 

' Compute the mean radius squared of the free surface at each node: 

CALL Rad (CoordS () , RadCurv2()) 



' Compute the partial wave height for the trial functions WH = M*Potential 
' where M=the mass distribution = ConstMass - (S/Smax) ^ExpMass 
CALL WaveHite (CoordS () , ModeS () , WavHiteO) 

' Compute the derivatives: R<WH)' and (R(WH)')'. By numerical experi- 
' mentation, it has found that the derivative (R(WH)')' must be formed 
9 slightly differently than R(WH)'. 

' Compute the (WH)' derivative 

CALL Deriv(l, CoordS (), WavHiteO, WavDerivl 0 ) 

' Now multiply by R to get R(WH') 

FOR I = 1 TO NoModes% 

FOR J = 1 TO NoNodes% 

WavDerivl (I, J) = WavDerivl (I, J) * CoordS (2, J) 

NEXT J 
NEXT I 

' Compute the (R(WH)')' derivative 

CALL Deriv(2, CoordS 0 , WavDerivl() f WavDeriv2()) 

' Compute and store in Al() the j-th part of the AO integrals. This 
9 is the part that is multiplied by WH-j and then integrated. 

CALL APart (CoordS 0 , RadCurv2 () , WavDeriv2 () , WavHiteO, A1 () ) 

r Compute the total Aij integrals 

CALL AIntegrals (FreqO , A1 () , CoordS () , WavHiteO, AO) 

' Compute the Bij integrals 

CALL BIntegrals (CoordS () , ModeS () , WavHiteO, B()) 

' Compute the eigenf requencies and eigenvectors. We must find the 
' eigenfrequencies by trial and error since A() is not 

' symmetric, and so most eigenvalue extraction routines will not work. 

FirstGuess = 0: Delta = .1 

CALL Root (FirstGuess, Delta, A(), B(), Freq2) 

9 The first Eigenf requency Freq2 has been computed, so now compute 
' the eigenvectors Vec() (or modal participation factors) 

CALL Vector (Freq2, A{), B(), Vec 0 ) 

' ========= END OF SEGMENT 3 END OF SEGMENT 3 ================== 

' =========== SEGMENT 4 SEGMENT 4 SEGMENT 4 =================== 

9 This segment computes the parameters of the mechanical model. 

' First, calculate the volume of liquid (non-dimensional) : LiqVolume 

CALL LiqVol (CoordS () , CoordW(), LiqVolume, LiquidCG, LastDZsurf, DeltaElev, 

SurfElev) 

' Calculate the non-dimensional slosh mass (ratio to liquid mass that 
' fills the tank completely) . 

CALL SloshM (CoordS () , CoordWO, ModeSO, ModeW(), FreqO, Vec(), Freq2, SMass, 
MFRatio) 

9 Compute the other mechanical model parameters 

CALL ParamModel (MFRatio, Freq2, SMass, SMass2, PendL, Spring, HO) 

' ============= END OF SEGMENT 4 END OF SEGMENT 4 ========== 

' ========== SEGMENT 5 SEGMENT 5 SEGMENT 5 ================= 

' This segment displays the results on the terminal screen, and also prints 
' the results out if the user desires. The printed results can include 
' the details of the wave height amplitude along the surface. 

CALL OutputScreen (Freq2, SMass2, PendL, HO, LiquidCG, Spring, Vec(), 
PrintCode%) 



IF PrintCode% = 0 THEN END ' user did not want printout 

CALL Printout (Freq2, FreqO, SMass2, PendL, HO, LiquidCG, Spring, Vec ( ) , 
PrintCode%) 

IF PrintCode% 2 THEN 9 user also wanted wave height data 

CALL Print0ut2 (CoordS () , WavHite ( ) , Vec(), Freq2, FreqO) 

END IF 

SCREEN 0 9 restore a clean screen 

END 





AINTEGRALS ★★*★*★***★*■*•*★***★*****★** 

' This routine computes each entry in the Aij matrix by numerical 
' integration, using linear interpolation. F{) is the ADINA eigenvalues, 

' APartialO is the j-th part of each integrand APart (), Coord(3,I) is the 
' "S" coordinate of each node point, WHO is the wave height (not including 
' the freq. term), and the answers are stored in Result () to be returned 
' to the calling program. 

SUB AIntegrals (F(), APartialO, Coord (), WHO, Result 0 ) 

SHARED NoNodes%, NoModes% 

FOR I - 1 TO NoModes% 

FOR J = 1 TO NoModes% 

Sum = 0 
AFreq = F ( J) 

FOR K = 1 TO NoNodes% 

SELECT CASE K ' First determine the right Delta-S 

CASE 1 ' first delta-S interval 

Delt = (Coord(3, 2) - Coord(3, 1)) / 2 
CASE NoNodes% ' last delta-S interval 

Delt = (Coord(3, NoNodes%) - Coord(3, (NoNodes% - 1))) / 2 
CASE ELSE ' interior delta-S interval 

Delt = (Coord (3, (K + 1) ) - Coord(3, (K - 1) ) ) / 2 
END SELECT 

Sum = Sum + APartial(J, K) * WH(I, K) * Delt 
NEXT K 

Result (I, J) = Sum * AFreq 
NEXT J 
NEXT I 

END SUB 





APART **************************** 

' This subroutine computes the j-th part of the Aij integral at each 
' node point, for later numerical integration. Coord (2, I) is the "R" 

' coordinate of each node and Coord (4,1) is R', Curvat2() is the square 
' of the mean curvature at each node, WD2() is the second derivative 
' computed by Deriv(2> at each node, WHO is the wave height (not 
' including freq) at each node, and the answer is stored in Result {) 

' to be sent back to the calling program. 

r 

SUB APart (Coord (), Curvat2() f WD2(), WH(), Result ()) 

SHARED NoModes%, NoNodes%, BondNo 

FOR I = 1 TO NoModes% 

FOR J = 1 TO NoNodes% 

R = Coord(2, J) : RPrime = Coord<4, J) : C2 = Curvat2 ( J) * Coord(2, J) 
IF J * 1 THEN 

Second = (C2 - BondNo * R * RPrime) * WH(I, J) 

ELSE 

Second = (C2 - 1 / R - BondNo * R * RPrime) * WH(I, J) 

END IF 

First = WD2 (I, J) 

Result (I, J) = First + Second 
NEXT J 
NEXT I 

END SUB 



r ***************************** BINTEGRALS ************************** 

' This subroutine computes the elements of the array Bij by numerical 
' integration, using linear interpolation. Coord (2, I) is the "R" 

' coordinate of each node, FuncO is the normalized potentials, WHO is 
' the wave height (not including the freq. term), and the answers are 
' stored in Result {) to be returned to the calling program. 

t 

SUB BIntegrals (Coord 0, FuncO, WHO, Result ()) 

SHARED NoNodes%, NoModes% 

FOR I - 1 TO NoModes% 

FOR J = 1 TO NoModes% 

Sum = 0 

FOR K * 1 TO NoNodes% 

SELECT CASE K 
CASE 1 

Delt * (Coord (3, 

CASE NoNodes% 

Delt = (Coord(3, 

CASE ELSE 

Delt = (Coord (3, 

END SELECT 

' Got to skip over the number of the node in FuncO 
Sum « Sum + Coord (2, K) * Func(J + 1, K) * WH(I, K) * Delt 
NEXT K 

Result (I, J) = Sum 
NEXT J 
NEXT I 

END SUB 


9 First compute the correct Delta-S 

' first Delta-S interval 
2) - Coord (3, 1) ) / 2 

9 last Delta-S interval 
NoNodes%) - Coord(3, (NoNodes% - 1))) / 2 
' interior Delta-S interval 
<K + 1)) - Coord (3, (K - 1) ) ) / 2 



r a****************************** DERIV **************************** 

' This subroutine computes the derivative of FuncDat ( ) and stores it in 
' FuncDerivO. Coord{3,I) is the "S" coordinate along the surface. 

' The derivatives are complicated because the spacing of the "S" points 
r is not uniform. A quadratic is fitted through each three points and 
' the derivative of the middle point is evaluated. Special formulas 
' are used for the first and last points. When Code%=2, the derivative 
' is the average of the linear (backwards) derivative and the quadratic derivative 

r 

SUB Deriv (Code%, CoordO, FuncDat (), FuncDerivO) 

SHARED NoModes%, NoNodes% 

IF Code% = 1 THEN 

AvgFactl = 0: AvgFact2 = 1 ' Discard the linear average 

ELSE 

AvgFactl = .5: AvgFact2 = .5 ' Weight the linear and quad derivs. 

equally 
END IF 


FOR I = 1 TO NoModes% 

FOR J « 1 TO NoNodes% 
SELECT CASE J 

CASE 1 

DF = Coord (3, 
DR = CoordO, 
CASE NoNodes% 

DF = Coord (3, 
DR = CoordO, 
CASE ELSE 

DF *= CoordO, 
DR = Coord(3, 
END SELECT 
SELECT CASE J 


' First, compute the backward and forward 
' differences of S 

2) - CoordO, 1) 

2 ) 

NoNodes%) - CoordO, (NoNodes% - 1)) 

(NoNodes% - 1) ) - CoordO, (NoNodes% - 2) ) 

(J + 1) ) - CoordO, J) 

J) - CoordO, (J - 1) ) 

Now compute the values of the function at the rear, 
middle, and forward point, and the appropriate 
differences 


CASE 1 

FR = -FuncDat (I, 2) 

FF = FuncDat (I, 2): FM = FuncDatd, 1) 

First = (FM - FR) / DR 

Second = FM * (DF - DR) / (DF * DR) 

Third = (FF * DR / DF - FR * DF / DR) / (DF + DR) 

FuncDerivO, 1) = First * AvgFactl + (Second + Third) * AvgFact2 
CASE NoNodes% 


FF = FuncDat(I, NoNodes%) 

FM « FuncDatd, (NoNodes% - 1)) 

FR = FuncDatd, (NoNodes% - 2)) 

First = (FF - FM) / DF 

Second = FM * (DF + DR) / (DF * DR) 

Third = (FF * (2 + DR / DF) + FR * DF / DR) / (DF + DR) 

FuncDeriv (I, NoNodes%) = First / 2 + (Third - Second) / 2 

CASE ELSE 

FF = FuncDatd, (J + 1)) 

FM = FuncDatd, J) 

FR = FuncDatd, (J - 1)) 

First = (FM - FR) / DR 

Second = FM * (DF - DR) / (DF * DR) 

Third = (FF * DR / DF - FR * DF / DR) / (DF + DR) 

FuncDeriv ( I, J) = First * AvgFactl + (Second + Third) * AvgFact2 
END SELECT 
NEXT J 
NEXT I 


END SUB 



r a************************** INPUT SCREEN **★*★*★***★*★*★*★*★***★*★**** 

' Sets up the terminal screen for input. This uses a QUICKPAK routine 
' Text Ini () that allows the arrows to be used to input data anywhere and 
' to back up through the data to correct any bad entries. If QUICKPAK 
' is not available, the data can be entered (after modifying the code) 

' but it will be more difficult to correct a bad entry without starting 
9 over from scratch for each set of input data. 

' Meaning of the parameters: 

' File%: 0 = data will be typed in from the terminal 

' 1 = data will be read in from the data disk files 

' Code% : 1 = input surface node # and R and S 

' 2 = input surface R f and Z ' 

' 3 = input surface R" and Z" 

' 4 = input surface potential values for each mode 

r 5 = input wall node # and R and Z 

' 6 = input wall potential values for each mode 

' Temp$ { ) : used to transfer data to Txt$ for display on the 

' screen, and to store input data temporarily after 

' it has been entered 

' TempNoO: used to accumulate the input data and return the 

' data to the main program for storage in the right 

' array 

' FileCode%: 1 = "open” disk files for data input 

' 2 = "append" data to disk files 

SUB Input Screen (File% / Code%, Temp$ <) , TempNoO) 

SHARED NoModes%, ModeNo%, NoNodes% 

CLS : COLOR 11, 0: FileCode% - 0 

' ============ SECTION 1 SECTION 1 SECTION 1 ================== 

' This section puts the right titles and instructions on the screen. 

SELECT CASE Code% 

CASE 1 TO 4 

LOCATE 1 , 20: 

LOCATE 1 , 29: 

LOCATE 1, 42: 

CASE 5 TO 6 

LOCATE 1, 20: 

LOCATE 1 , 29: 

LOCATE 1 , 39: 

END SELECT 

LOCATE 2, 10 
SELECT CASE Code% 

CASE 1 

PRINT "Enter: Node #, and R, S coordinates of node, starting at R=0. 

CASE 2 

PRINT "Enter: dR/dS and dZ/dS derivatives for each node." 

CASE 3 

PRINT "Enter: d2R/dS2 AND d2Z/dS2 derivatives of each mode.” 

CASE 4 

PRINT "Enter: Potential values for mode number " 

LOCATE 2, 49: COLOR 14, 0: PRINT ModeNo% 

CASE 5 

PRINT "Enter: Node #, and R, Z coordinates of node, starting at R=0" 

CASE 6 

PRINT "Enter: Potential values for mode number " 

LOCATE 2, 49: COLOR 14, 0: PRINT ModeNo% 

END SELECT 


COLOR 

15, 

0: 

PRINT 

PRINT 

"DATA FOR " 
"FREE SURFACE " 

COLOR 

11, 

0: 

PRINT 

"NODES" 

COLOR 

15, 

0: 

PRINT 

PRINT 

"DATA FOR " 
"TANK WALL " 

COLOR 

11, 

0: 

PRINT 

"NODES" 



instructions for typing input data 


IF File% “ 0 THEN 
COLOR 15, 0 


LOCATE 4, 

8: 




PRINT 

”Use the ” 

COLOR 14, 

0: 

LOCATE 

4, 

16 

PRINT 

"ENTER ” 

COLOR 15, 

0: 

LOCATE 

4, 

22 

PRINT 

"or ” 

COLOR 14, 

0: 

LOCATE 

4, 

25 

PRINT 

"DOWN ARROW " 

COLOR 15, 

0: 

LOCATE 

4, 

37 

PRINT 

"key after entering each value.” 

LOCATE 5, 

8: 




PRINT 

"Use the ” 

LOCATE 5, 

16 

: COLOR 

14, 

0 

PRINT 

"UP ARROW " 

COLOR 15, 

0: 

LOCATE 

5, 

25 

PRINT 

"key to move back through the input 

data . " 






ELSE 





' instructions for reading input from files 

COLOR 15, 

0 






LOCATE 4, 

8: 




PRINT 

"Use the " 

COLOR 14, 

0: 

LOCATE 

4, 

16 

PRINT 

"ENTER " 

COLOR 15, 

0: 

LOCATE 

4, 

22 

PRINT 

"or ” 

COLOR 14, 

0: 

LOCATE 

4, 

25 

PRINT 

"DOWN ARROW " 

COLOR 15, 

0: 

LOCATE 

4, 

37 

PRINT 

"key to enter each value from the file." 

LOCATE 5, 

8: 




PRINT 

"Use the " 

LOCATE 5, 

16 

: COLOR 

14, 

0: 

PRINT 

"UP ARROW " 

COLOR 15, 

0: 

LOCATE 

5, 

25: 

PRINT 

"key to move back up through the data." 


END IF 


SELECT CASE Code% 


CASE 1 ' Print titles over input columns for R,S input data 



LOCATE 

7, 

5 : 

PRINT " 

Node 

# 


R 

s 



Node 

# 


R S" 














CASE 2 



' Print 

titles over 

input 

columns 

for R',Z' 

input 

data 



- 

LOCATE 

7, 

5: 

PRINT " 

Node 

# 


R' (S) 

Z' (S) 


Node 

# 


R' (S) Z' 

(S) " 













CASE 3 



' Print 

titles over 

input 

columns 

for R" , Z 

" input 

data 




LOCATE 

7, 

5: 

PRINT " 

Node 

# 


R' ' (S) 

Z ' ' (S) 


Node 

# 


R' ' (S) Z' 

' (S) 

ii 












CASE 4 



9 Print 

titles over 

input 

columns 

for nodal 

values 

input 




LOCATE 

7, 

5: 

PRINT " 

Node 

# 



Potential 

Node 

# 


Potential" 













„ _ . 

CASE 5 














LOCATE 

7, 

5: 

PRINT " 

Node 

# 


R 

Z 



Node 

# 


R Z 

ri 













CASE 6 



' Print 

titles over 

input 

columns 

for nodal 

values 

input 



— 

LOCATE 

7, 

5: 

PRINT " 

Node 

# 



Potential 

Node 

# 


Potential” 

END SELECT 

' ========= end OF SECTION 1 END OF SECTION 1 ============= 

' ============== SECTION 2 SECTION 2 SECTION 2 ============= 

' This section gets the screen display in shape to continue with data 
' input . 

' Print out the numbers from 1 to NoNodes% on the screen in two 
' columns, 1 to 12 and 13 to NoNodes%, to help organize the input. 

FOR N s 1 TO NoNodes% 

IF N <= 12 THEN 

LOCATE 7 + N, 3: PRINT N 
ELSE 

LOCATE 7 + N - 12, 42: PRINT N 
END IF 
NEXT N 



' When entering the modal values, the surface or wall node numbers are 
' also displayed in the columns, but we have to retrieve them from Temp$ 

' and store them in TempNoO so they can be displayed. The node numbers 
' are already available when entering the coordinate data. 

IF Code% =* 4 OR Code% = 6 THEN 
FOR M * 1 TO NoNodes% 

TempNod, M) = VAL (Temp$ (1, M) ) 

NEXT M 
END IF 

' Now, for the wall input, change Code% so that the same routines can 
' be used as for the surface nodes; also change FileCode% so that data is 
' appended to the existing disk files rather than overwriting them. 

IF Code% - 5 THEN Code% - 1: FileCode% * 1 

IF Code% - 6 THEN Code% = 4: FileCode% * 1 

9 Now print the node numbers on the screen: 

SELECT CASE Code% 

CASE 2 TO 4 

FOR N - 1 TO NoNodes% 

IF N <= 12 THEN 

LOCATE 7 + N, 10: PRINT TempNo(l, N) 

ELSE 

LOCATE 7+N-12, 50: PRINT TempNod, N) 

END IF 
NEXT N 
END SELECT 


' ============= end OF SECTION 2 END OF SECTION 2 ========== 

' =========== SECTION 3 SECTION 3 SECTION 3 ================== 

' This long segment gets the screen input and checks to see if it is 
9 okay. Most of the code just allows the cursor to move up and down 

9 on the screen to correct any bad entry, without having to enter all 

9 the data over again. The parameter ExC% indicates whether the piece 
' of input is okay (ExC%=0) so the cursor can be moved to the next line, 

' or that the piece of data is wrong (ExC%<>0) and will be re-entered. 

9 ExC% is set by how the user ends the entry (RETURN, etc., or UP ARROW, etc) 

9 First, decide what column to put the cursor into initially. 

N = 1: NN = N 

IF Code% = 1 THEN Ml = 10 ' Put cursor in "Node No" col 

IF Code% = 2 OR Code* = 3 THEN Ml = 20 9 Put cursor in 2nd col 

IF Code% = 4 THEN Ml = 30 9 Put cursor in 3rd col 

' When typing in the modal values, we need to clear out the coordinate data 
' from Temp$ { ) so that it is not displayed on the screen incorrectly. 

IF Code% = 4 AND File% = 0 THEN 
FOR 1=1 TO NoNodes% 

Temp$ (ModeNo% + 1, I) - "" 

NEXT I 
END IF 

Insurf : ' we keep coming back here for the next piece of data 

COLOR 7 , 0 
SELECT CASE Ml 

' Txt$ is the data that will be displayed on the screen initially 
CASE 10, 50 

Txt$ - Temp$(l, N) 9 node numbers 

CASE 20, 60 

Txt$ = Temp$(2 * Code%, N) ' R, R' or R" (surface) or R (wall) 

CASE 30, 70 

IF Code* <> 4 THEN 

Txt $ = Temp$(l + 2 * Code*, N) f S, Z' or Z" (surface) or Z (wall) 
ELSE 



' modal potential values at nodes 


Txt$ = Temp$ (ModeNo% + 1, N) 

END IF 
END SELECT 

LOCATE 7 + NN, Ml 
' Now get the piece of data; it is returned in Txt$ 

CALL Text Ini (Txt$, 8, 0, 0, ExC%, 7) 

IF ExC% - 0 THEN 

' data is okay, so store it, move the cursor, and get the next data piece 
SELECT CASE Ml 
CASE 10, 50 

TempS (1, N) = Txt$ 

TempNo(l, N) = VAL(Txt$) ’ change the text input to a number 

Ml = Ml + 10 
GOTO Insurf 
CASE 20, 60 

Temp$ (2 * Code%, N) = Txt$ 

TempNo (2 * Code%, N) = VAL(Txt$) 

Ml = Ml + 10 
GOTO Insurf 
CASE 30 

IF Code% <> 4 THEN 

TempS (1 + 2 * Code % , N) = Txt$ 

TempNo (1 + 2 * Code%, N) = VAL(Txt$) 

ELSE 

Temp$ (ModeNo% + 1, N) = Txt$ 

TempNo (ModeNo% + 1, N) = VAL(Txt$) 

END IF 

IF N < 12 THEN 

N = N + 1; NN-N 
Ml = 10 

IF Code% = 2 OR Code% = 3 THEN Ml = 20 
IF Code% = 4 THEN Ml = 30 
ELSE 

N = N + 1: NN = N - 12 
Ml = 50 

IF Code% = 2 OR Code% = 3 THEN Ml = 60 
IF Code% = 4 THEN Ml = 70 
END IF 
GOTO Insurf 
CASE 70 

IF Code% <> 4 THEN 

Temp$ (1 + 2 * Code%, N) = Txt$ 

TempNo (1 + 2 * Code%, N) = VAL(Txt$) 

ELSE 

TempS (ModeNo% + 1, N) = Txt$ 

TempNo (ModeNo% + 1, N) = VAL(Txt$) 

END IF 

IF N < NoNodes% THEN 

N = N + 1: NN = N - 12 
Ml = 50 

IF Code% = 2 OR Code% = 3 THEN Ml = 60 
IF Code% = 4 THEN Ml = 70 
GOTO Insurf 
ELSE 

LOCATE 23, 10: COLOR 15, 0 

INPUT "Is all your data correct (Y, N) ? ", Responses 
Responses = UCASES (Responses) 

IF Responses = "Y" THEN 
GOTO Finish 

ELSE ' start this screen all over again 

N = NoNodes% 

NN = NoNodes% - 12 



put cursor on last item of data 


Ml = 70 
GOTO Insurf 
END IF 
END IF 
END SELECT 

ELSE 

' piece of data is bad, so don't store it and back up the cursor. This 
' is complicated because we have to jump to a new column when the cursor 
' has backed all the way to the top of the previous one, and we also 
' have to make sure that the cursor does not leave the data input area. 
SELECT CASE Ml 
CASE 30, 70 

IF Code% <> 4 THEN 
Ml = Ml - 10 
GOTO Insurf 
ELSE 

IF Ml = 70 THEN 
IF N * 13 THEN 
Ml = 30 
N = N - 1 
NN * N 
GOTO Insurf 
ELSE 

Ml - 70 
N = N - 1 
NN = N - 1 2 
GOTO Insurf 
END IF 
ELSE 

N = N - 1 
IF N < 1 THEN 
CALL Chime (8) 

N = 1 
NN = N 
GOTO Insurf 
ELSE 

NN * N 
GOTO Insurf 
END IF 
END IF 
END IF 
CASE 20 

IF Code% = 1 THEN 
Ml « Ml - 10 
GOTO Insurf 
ELSE 

N “ N - 1 
IF N < 1 THEN 
CALL Chime (8) 

N = 1 
Ml = 20 
GOTO Insurf 
ELSE 

Ml = 30 

NN = N 

GOTO Insurf 
END IF 
END IF 
CASE 60 

IF Code% = 1 THEN 

Ml = Ml - 10 

GOTO Insurf 



ELSE 

IF N = 13 THEN 
Ml = 30 
N = N - 1 
NN = N 
GOTO Insurf 
ELSE 

Ml = 70 
N = N - 1 
NN - N - 12 
GOTO Insurf 
END IF 
END IF 
CASE 10 

N = N - 1 
IF N < 1 THEN 
CALL Chime (8) 

N = 1 
Ml = 10 
GOTO Insurf 
ELSE 

Ml = 30 
NN * N 
GOTO Insurf 
END IF 
CASE 50 

IF N = 13 THEN 
Ml = 30 
N = N - 1 
NN = N 
GOTO Insurf 
ELSE 

Ml = 70 
N - N - 1 
NN * N - 12 
GOTO Insurf 
END IF 
END SELECT 
END IF 

' ============= END OF SECTION 3 END OF SECTION 3 ============ 

' ========== SECTION 4 SECTION 4 SECTION 4 ============== 

' This section writes the good data to the disk data files. 

Finish : 

' Write the surface coordinate data after all three screens of data 
' have been entered: 

IF Code% * 3 THEN 

OPEN "COORD.DAT" FOR OUTPUT AS #1 ' This wipes out all existing data 

FOR I * 1 TO 7 

FOR J - 1 TO NoNodes% 

WRITE #1, TempNo ( I , J) 

NEXT J 
NEXT I 
CLOSE #1 
END IF 

' Append the wall coordinate data after it has been completely entered: 

IF Code% = 1 AND FileCode% = 1 THEN 
OPEN "COORD.DAT” FOR APPEND AS #1 
FOR I = 1 TO 3 

FOR J * 1 TO NoNodes% 

WRITE #1, TempNo (I, J) 

NEXT J 



NEXT I 
CLOSE #1 
END IF 

' Write the nodal data after all the modes have been entered: 

IF Code% = 4 AND ModeNo% « NoModes% THEN 

' check to see if the data is for the surface nodes: 

IF FileCode% - 0 THEN 

OPEN "MODE.DAT" FOR OUTPUT AS #1 ' This wipes out existing data 

FOR 1=1 TO NoModes% 

FOR J - 1 TO NoNodes% 

WRITE #1, TempNod + 1, J) ' Skip the node #'s 

NEXT J 
NEXT I 
CLOSE #1 
END IF 

' check to see if the data is for the wall nodes: 

IF FileCode% = 1 THEN 

OPEN "MODE.DAT" FOR APPEND AS #1 
FOR I = 1 TO NoModes% 

FOR J = 1 TO NoNodes% 

WRITE #1, TempNo(I + 1, J) ' Skip the node #'s 

NEXT J 
NEXT I 
CLOSE #1 
END IF 
END IF 

END SUB 





LIQVOL **************************** 

' This subroutine computes the non-dimensional volume occupied by the 
' liquid and the non-dimensional location of the center of mass. 

' The volume = pi * RT A 3 * K, where RT is the non-dimensional 
' radius to the tank wall from the origin of the coordinate system. 

' RT must be equal to ONE if the input is given correctly. This routine 
' computes K, assuming that RT=1. The dimensional volume is the non- 
' non-dimensional volume multiplied by the radius RO to the tank wall; 

' or Volume = pi * R0 A 3*K. The center of mass is referenced to the 
9 bottom of the tank. 

r 

SUB LiqVol (Surf () , Wall(), Ansi, Ans2, DZ, DeltaElev, SurfElev) 

SHARED NoNodes% 

DIM ZSurf (25) 

' We first have to compute the Z coordinates of the free surface nodes 
' since they were not asked for in the input. The Z's are computed 
' by integrating the values of Z' . 

ZSurf (1) = 0 

FOR I = 2 TO NoNodes% 

AvgSlope - (Surf (5, (I - 1) ) + Surf (5, I) ) / 2 
DeltaZ = AvgSlope * (Surf (3, I) - Surf (3, (I - 1))) 

ZSurf (I) = ZSurf (I - 1) + DeltaZ 
NEXT I 

DZ = ZSurf (NoNodes% ) - ZSurf ( (NoNodes% - 1)) 

' Estimate the elevation of the free surface center from the tank bottom 
' First, compute distance tank bottom is below Z = 0: 

DeltaElev - -Wall (3, 1) 

9 Then, the free surface elevation is: 

ZSurf (1) = DeltaElev + (Wall (3, NoNodes%) - ZSurf (25)) 

' Save ZSurf ( 1 ) 

SurfElev = ZSurf (1) 

' Correct the other free surface elevations: 

FOR I = 2 TO NoNodes% 

ZSurf (I) = ZSurf ( I ) + ZSurf(l) 

NEXT I 

9 Now can do the numerical integration. Differential volume = 

' (Rad to wall) A 2 * DeltaZ of wall - 

' (Rad to liq surface) A 2 * DeltaZ of surface 

Ansi = 0: Ans2 = 0 

An si Temp = 0: An s 2 Temp = 0 

An s 3 Temp - 0: An s 4 Temp = 0 

FOR I = 1 TO NoNodes% 

9 First set up the right DeltaZ' s 
SELECT CASE I 

CASE 1 ' first Delta-Z interval 

DZSurf = (ZSurf (2) - ZSurf (1)) / 2 

DZWall = (Wall (3, 2) - Wall (3, 1) ) / 2 

CASE NoNodes% ' last Delta-Z interval 

DZSurf = (ZSurf (NoNodes%) - ZSurf (NoNodes% -1)) / 2 
DZWall = (Wall (3, NoNodes%) - Wall (3, (NoNodes% -1))) / 2 
CASE ELSE ' interior Delta-Z interval 

DZSurf * (ZSurf (I + 1) - ZSurf (I -1)) / 2 
DZWall = (Wall (3, (I + 1) ) - Wall (3, (I - 1) ) ) / 2 
END SELECT 



DeltaVoll - (Wall (2, I) A 2) * DZWall 
DeltaVol2 = (Surf (2, I)) - 2 * DZSurf 
An s 1 Temp = AnslTemp + DeltaVoll 
Ans2Temp - Ans2Temp + DeltaVol2 

Ans3Temp * Ans3Temp + (Wall (3, I) + DeltaElev) * DeltaVoll 
Ans4Temp = Ans4Temp + ZSurf(I) * DeltaVol2 
NEXT I 

Ansi “ AnslTemp - Ans2Temp ' liquid volume 

Ans2 * (Ans3Temp - Ans4Temp) / Ansi ' liquid c.g. 

END SUB 



r ************************** OUTPUTSCREEN ************************** 

' This subroutine displays the input data and the computed results on 
' the screen. It also asks if the user wants printed copies. 

f 

SUB OutputScreen (F2, SM2, PL, HO, CG, SC, BI(), PCode%) 

SHARED FillLevel, BondNo, Ang, ConstMass, ExpMass, NoModes%, NoNodes% 
SHARED LiqVolume 

Pi = 4 * ATN(l) 

CLS 

' ================= Display the input ============================ 

COLOR 11: LOCATE 1, 35: PRINT "INPUT DATA" 

COLOR 15: LOCATE 2, 5: PRINT "Bond Number Fill Level % Cont . Ang" 

LOCATE 2, 44: PRINT "Surface Mass Distribution Function" 

CALL Box (3, 7, 5, 13, 1, 14, -1) ' draw some boxes for the output 


CALL 

Box (3, 

20, 

5, 

27, 

1, 

14, 

-1) 

CALL 

Box (3, 

33, 

5, 

40, 


14, 

-1) 

CALL 

Box (3, 

48, 

5, 

74, 

1, 

14, 

-1) 


LOCATE 

4, 

8: 

PRINT 

USING 

BondNo 

' format and print output 

LOCATE 

4, 

22: 

PRINT 

USING 

FillLevel 

LOCATE 

4, 

35: 

PRINT 

USING 

Ang 


LOCATE 

4, 

50: 

PRINT 

USING 

"## . ###" ; ConstMass 


LOCATE 

4, 

56: 

PRINT 

USING 

" - (S/Smax) 

ExpMa s s 


================ end of input display ======== 


' ======================= display the computed results ============== 

COLOR 11: LOCATE 6, 22: PRINT "PARTICIPATION FACTORS OF TRIAL MODES" 

FOR I = 1 TO 3 

Edgel% = 10 + 25 * (I - 1) 

Edge2% = Edgel% + 10 

CALL Box (7, ( Edge 1 % ) , 9, (Edge2%), 1, 14, -1) 

NEXT I 

COLOR 15: LOCATE 8, 4 


PRINT "B (1) =": 
LOCATE 8 , 29 

LOCATE 

8, 

12 : 

PRINT 

"1.00000" 

PRINT "B(2) =": 
LOCATE 8, 54 

LOCATE 

8, 

37: 

PRINT 

USING "#.#####"; BI (2 ) 

PRINT "B(3) =": 
COLOR 11 

LOCATE 

8, 

62: 

PRINT 

USING "#.#####"; BI (3) 


LOCATE 10, 15: PRINT "NON-DIMENSIONAL PARAMETERS OF THE PENDULUM MODEL" 

COLOR 15 

LOCATE 11, 2: PRINT " (Liq Vol)/("; CHR$(227); "* Ro^) " 

LOCATE 11, 27: PRINT "Slosh M/Liquid M" 

LOCATE 11, 49: PRINT "Pend. L/Ro" 

LOCATE 11, 68: PRINT "Ho/Ro" 

CALL Box ( 12 , 5, 14, 15, 1, 14, -1) 

LOCATE 13, 7: PRINT USING "##.####"; LiqVolume 

CALL Box (12, 29, 14, 39, 1, 14, -1) 

LOCATE 13, 31: PRINT USING "#.####"; SM2 

CALL Box (12, 48, 14, 59, 1, 14, -1) 

LOCATE 13, 50: PRINT USING "##.####"; PL 

CALL Box (12, 66, 14, 75, 1, 14, -1) 

LOCATE 13, 67: PRINT USING "###.###"; HO 

LOCATE 15, 5: PRINT "CG Loc./Ro" 

LOCATE 15, 28: PRINT "Spring/ ("; CHR$(229); " * Ro / '2) " 

LOCATE 15, 50: PRINT "Freq / '2/("; "(1+Bo)"; CHR$(229); "/(d * Ro A 3))" 

CALL Box (16, 5, 18, 16, 1, 14, -1) 

CALL Box (16, 26, 18, 45, 1, 14, -1) 

CALL Box (16, 58, 18, 71, 1, 14, -1) 



LOCATE 17, 7: PRINT USING "##.####"; CG 

LOCATE 17, 32: PRINT USING ”##.####"; SC 

LOCATE 17, 62: PRINT USING ”##.####"; F2 


end of results display ===================== 


LOCATE 21, 5 

PRINT "Do you want the above data printed out (Y or N) ?" 

LOCATE 21, 55, 1, 5, 7: INPUT Responses 

Responses ** UCASES (Responses) 

IF Responses - "N" THEN 
PCode% - 0: EXIT SUB 

ELSE 

PCode% = 1 
END IF 

LOCATE 23, 5 

PRINT "Do you also want the surface wave data printed (Y or N) ?" 
LOCATE 23, 63: INPUT "", Responses 
Responses = UCASES (Responses) 

IF Responses = "N" THEN PCode% = 1 ELSE PCode% = 2 
END SUB 





★★★★★★★★★★★★★★★★A************ PAD ZERO 
This subroutine turns a number into a text string for printing. The string 
can be specified to have a given number of decimal places "NoDecs" and 
a given number of digits before the decimal "NoBef oreDec" If the 
number is a floating point number with a negative exponent and the 
value of the exponent is greater than the NoDecs, then the string 
is set equal to zero; otherwise the number is turned into a string that 
looks like a decimal number. At this time, nothing is done to 
floating point numbers with positive exponents; the string is just the 
number. 


FUNCTION PadZero$ (A, NoDecs, NoBeforeDec) 

A$ = STR$ (A) 

A$ = LTRIM$ (RTRIM$ (A$ ) ) ' Get rid of all blank spaces 

C$ = MID$ (A$, 1, 1) 

' In case the number is negative, keep the same number of digits before 
' the decimal point. 

IF C$ * THEN NoBeforeDec = NoBeforeDec + 1 


9 find if it is a floating point number 

FOR 1 = 2 TO LEN(A$) ' skip any possible signs for negative #'s 

B$ = MID$ ( A$ , I, 1) 

IF B$ = THEN Flag = 1: NF = 1 + 1 9 negative exponent 

IF B$ = "+” THEN Flag = 2: NF = 1 + 1 9 positive exponent 

NEXT I 


SELECT CASE Flag 
CASE 1 

Numl = VAL (MID$ ( A$ , NF, 1) ) : 
IF Numl > 0 OR Num2 > NoDecs 
A$ = ”0." 

FOR I = 1 TO NoDecs 

A$ = A$ + "0” 

NEXT I 
ELSE 

IF C$ = THEN 

B$ = MID$(A$, 2, 1) 

FOR I = 4 TO (LEN (A$ ) 
B$ = B$ + MID$(A$, I, 
NEXT I 

FOR I = 1 TO (Num2 - 1 
B$ = ”0” + B$ 

NEXT I 

A$ = C$ + " 0 . ” + B$ 
ELSE 

B$ = MID$(A$, 1, 1) 

FOR I = 3 TO ( LEN (A$ ) 
B$ = B$ + MID$ ( A$ , 
NEXT I 

FOR I = 1 TO (Num2 - 1 
B$ = ”0” + B$ 

NEXT I 

A$ = "0." + B$ 

END IF 
END IF 

CASE ELSE ' The n 

END SELECT 


The number has a negative exponent 

Num2 = VAL (MID$ (A$, (NF + 1), 1) ) 
THEN 9 the number is too small 

9 make it equal to 0.000... 


turn the number into a decimal number 


4) 


1 ) 

) 


get the first digit 
skip - sign, 1st digit, and 
the decimal point, and don't 
count the M E-XX" at end 


- 4) 
I, 1 ) 

) 


get the first digt 

skip 1st digit and decimal point 

and don't count the "E-XX M 


umber has no exponent or a postive exponent 


' if it is a negative number, make it start with "-0." 
IF C$ = AND MID$(A$, 2, 1) = " . " THEN 

B$ = "" 

FOR I = 3 TO LEN (A$ ) 

B$ = B$ + MID$ ( A$ , I, 1) 



NEXT I 

A$ = "-0." + B$ 

END IF 

' Find out where the decimal point is 

DecLoc = 0 

FOR I = 1 TO LEN (A$ ) 

B$ « MID$ (A$, I, 1) 

IF B$ - n . " THEN DecLoc * I 
NEXT I 

' if the first character is ” . " then add a zero before it 
IF DecLoc * 1 THEN A$ - "0” + A$ : DecLoc = 2 
IF C$ * AND MID$<A$, 2, 1) - " . " THEN 

B$ = "" 

FOR I = 3 TO LEN (A$) 

B$ * B$ + MID$ ( A $ f I, 1) 

NEXT I 

A$ - "-0." + B$ 

END IF 

' if there is no " ♦ " then add one at the end 

IF DecLoc = 0 AND NoDecs > 0 THEN A$ = A$ + " . " : DecLoc * LEN (A$) 

' if the number has too many decimal places, then cut the extra ones 
IF ( (LEN (A$) - DecLoc) >* NoDecs) THEN 
B$ = 

FOR I = 1 TO (DecLoc + NoDecs) 

B$ = B$ + MID$ ( A$ r I, 1) 

NEXT I 
A$ = B$ 

END IF 

' if there are not enough ”0” at the end, add enough of them, but 
' first, find out how many "0" to add 
NoPad * NoDecs - (LEN(A$) - DecLoc) 

' now add the "0"s 
FOR I ■ 1 TO NoPad 
A$ = A$ + ”0” 

NEXT I 

' now see if there are enough spaces before the decimal point 
IF NoBeforeDec > (DecLoc - 1) THEN ' if not, add blank spaces 
NoPad = NoBeforeDec - (DecLoc - 1) 

FOR I - 1 TO NoPad 
A$ * " " + A$ 

NEXT I 
END IF 

9 All done, so send the string back 
PadZero$ = A$ 


END FUNCTION 



' *************************** PARAMMODEL ************************ 

' This subroutine computes the non-dimensional parameters of the 
' pendulum model: SM2 = Slosh Mass / Liquid Mass in tank; 

' PL = Pendulum Length / RO 

' SC = Spring Constant / (surface tension * RCTS) 

' HO = HingePointLocat ion / RO 

' Note: F2 = SlsFreq (i.e, non-dimensional natural frequency * 2) 

' SM = SloshMass / LiquidMass in Full Tank 

' MOF = Moment / Force (non-dimensional) 

r 

SUB ParamModel (MOF, F2, SM, SM2, PL, SC, HO) 

SHARED LiqVolume, FillLevel, DeltaElev, SurfElev 
SHARED BondNo 

Pi = 4 * ATN(l) 

' Compute the slosh mass ratio to the actual liquid mass 
SM2 = SM / (FillLevel / 100) 

' Compute the pendulum length ratio to tank radius 
PL = 1 / F2 

' Compute the spring constant ratio to (surface tension * RO A 2) 

SC = SM2 * Pi * LiqVolume / F2 

' Compute the hinge point location ratio to tank radius 

HO = MOF + SC / ((1 + BondNo) * (SM2 * Pi * LiqVolume) * PL * F2) 

HO = HO + SurfElev ' reference to tank bottom 

END SUB 



r ************************* POTNORM ******************************* 

9 This subroutines finds the largest value of the trial potentials for 
9 each mode and uses it to normalize all the nodal values for that mode. 

t 

SUB PotNorm (SurfPO, WallPO) 

SHARED NoModes%, NoNodes% 

FOR 1 = 2 TO (NoModes% + 1) r Skip over the node ID number 

PotMax = ABS(SurfP(I, 1)) ' First guess for maximum value 

FOR J « 2 TO NoNodes% 

IF ABS (Surf P ( I , J) ) > PotMax THEN PotMax = ABS(SurfP<I, J) ) 

NEXT J 

FOR J = 1 TO NoNodes% 

SurfP(I, J) = SurfP(I, J) / PotMax ' This is the normalization 

WallP(I, J) = WallP (I, J) / PotMax 
NEXT J 
NEXT I 

END SUB 



r ***************************** PRINTOUT **************************** 

' This printout routine prints out the input data, the trial results, 

' and the mechanical model paramters. It uses the function "PadZero$ M 
' to make strings with a specified number of decimal places and leading 
' digits out of the numerical data for easier control of the printing. 

9 

SUB Printout (F2, F(), SM2, PL, HO, CG, SC, BI(), PC%) 

SHARED FillLevel, BondNo, Ang, ConstMass, ExpMass, NoModes%, NoNodes% 
SHARED LiqVolume 

DIM F$(3) 

CLS 


' ================= Display the input ============================ 

LPRINT SPC (27) ; "===== INPUT DATA =====" 

A$ = "Bond Number Fill Level (%) Cont . Ang" 

B$ - "Surface Mass Distribution Function" 

LPRINT SPC ( 4 ) ; A$ ; SPC ( 4 ) ; B$ 

B$ = PadZero$ (BondNo, 1, 2) '1 decimal place, two spaces or digits 

F$ = PadZero$ (FillLevel, 1, 2) ' (cont'd) before the decimal place 

A$ = PadZero$ (Ang, 1, 2) 

C$ = PadZero$ (ConstMass, 3, 2) 

E$ = PadZero$ (ExpMass, 1, 2) 

LPRINT SPC (7 ) ; B$; SPC(IO); F$; SPC(9); A$; SPC(12); C$; " - (S/Smax)''"; E$ 
LPRINT 

LPRINT SPC (20) ; 

LPRINT "Trial Freq 1"; SPC(4); "Trial Freq 2"; SPC(4); "TrialFreq 3" 

FOR I ■ 1 TO 3 

F$ ( I ) = PadZero$( (F(I) ) , 5, 2) 

NEXT I 

LPRINT SPC (13) ; 

FOR I = 1 TO 3 
LPRINT SPC (8); F$(I); 

NEXT I 

' ====================== end of input display ==================== 

' ======================= display the results ===================== 

LPRINT : LPRINT SPC (27); "======= RESULTS =======" 

LPRINT SPC (15); "***** PARTICIPATION FACTORS OF TRIAL MODES *****" 

FOR I = 1 TO 3 

F$ ( I ) = PadZeroS ( (BI (I) ) , 5, 2) 

NEXT I 

LPRINT SPC (10) ; 

FOR I = 1 TO 3 

1$ = LTRIM$ (RTRIM$ (STR$ (I) ) ) 

LPRINT SPC (4) ; "B("; 1$; ") = "; F$(I); 

NEXT I 

LPRINT : LPRINT 
LPRINT SPC ( 10 ) ; 

LPRINT "***** NON-DIMENSIONAL PARAMETERS OF THE PENDULUM MODEL *****" 

L$ = PadZeroS (LiqVolume, 4, 2) 

S$ - PadZeroS (SM2, 4, 1) 

P$ = PadZero$ (PL, 4, 2) 

H$ = PadZero$ (HO, 3, 3) 

LPRINT SPC (10) ; 

LPRINT "Liq Vol/"; CHR$(227); "Ro^3"; SPC(3); 

LPRINT "Slosh Mass/Liq Mass"; SPC(3); 

LPRINT "Pend. L/Ro"; SPC(7); 

LPRINT "Ho/Ro" 

LPRINT SPC (12); L$; SPC (13); S$; SPC (11); P$; SPC(6); 


H$ 



C$ = PadZero$(CG, 4, 2) 

S$ - PadZero$ (SC, 4, 2) 

F$ = PadZero$(F2, 4, 2) 

LPRINT : LPRINT SPC(IO); 

LPRINT "Liq. CG/Ro"; SPC(8); 

LPRINT "Spring/"; CHR$(229); "Ro A 2) n ; SPC(8); 

LPRINT "Freq A 2/["; "(1+Bo)"; CHR$(229); "/(dRo A 3)]" 

LPRINT SPC(ll); C$, SPC(2); S$; SPC{20); F$ 

' ================= end of results display ===================== 

IF PC% = 1 THEN LPRINT CHR$(12) ELSE LPRINT : LPRINT 

END SUB 



******************************* PRINT0UT2 

This subroutine prints out the wave height data 




SUB Print0ut2 (S(), WH(), BI(), F2, F()) 

SHARED NoNodes%, NoModes%, BondNo 
DIM W ( 4 ) , W$ (4) 

LPRINT SPC (22 ) ; "| ******* WAVE HEIGHTS ON SURFACE ******* |" 
LPRINT SPC (23) ; "ADINA"; SPC(7); "ADINA" ; SPC(7); "ADINA"; 

LPRINT SPC (7); "Slosh" 

LPRINT SPC (10) ; "S Coord."; SPC(4); "Trial 1"; SPC(5); "Trial 2"; 
LPRINT SPC (5) ; "Trial 3"; SPC(6); "Mode 1" 

LPRINT SPC (9) ; " "; SPC (5) ; " SPC (5) ; ” "; 

LPRINT SPC (5) ; " "; SPC (5) ; " " 


FOR I = 1 TO NoNodes% 

LPRINT SPC (11) ; 

S$ = PadZero$ (S (3, I), 4, 1) 
W(l) = WH ( 1 , I) 

W (2) = 0 
W (3) = 0 


= WH (2 , I) 


WH (2, 
WH (3, 


I) 

I) 


BI(J) * F ( J) 
SQR (F ( J) ) * 


(1 


SELECT CASE NoModes% 

CASE 2 
W(2) 

CASE 3 
W (2 ) = 

W (3) = 

CASE ELSE 
END SELECT 
W(4) = 0 
FOR J = 1 TO 3 
W (4) = W { 4 ) 

W(J) = W(J) 

NEXT J 

W(4) = W ( 4 ) * (1 + BondNo) / SQR(F2) 
FOR J = 1 TO 4 

W$(J) = PadZeroS (W(J) 

NEXT J 
W = 5 
5 
5 
5 

= LTRIM$ (W$ (1) ) 

= LTRIM$ (W$ (2) ) 

= LTRIM$ (W$ (3) ) 

= LTRIM$ (W$ (4) ) 

S$ ; SPC (X) ; 


W(J) 

+ BondNo) 


4, 2) 


X = 

y = 

z = 

B$ 

B$ 

B$ 

B$ 
LPRINT 
NEXT I 


IF 

MID$ (B$, 

1, 

1) = then 

X = 4 

IF 

MID$ (B$ , 

1, 

1) = THEN 

Y = 4 

IF 

MID$ (B$ , 

1, 

1) = THEN 

Z = 4 

IF 

MID$ (B$ , 

1, 

1) = "->> THEN 

W — 4 


W$ ( 1 ) ; SPC(Y); W${2); SPC(Z); W$(3); 


SPC (W) ; 


LPRINT CHR$(12) ' feed paper out 

END SUB 


W$ (4) 





p fip *********★★*★***★*★*★*★*★*★*★ 

' This subroutine computes the expression (1/R1^2 + 1/R2 /V 2) where R1 
' and R2 are the mean radii of curvature of the free surface. R1 = 

' Z (S) ”R(S) ' -Z (S) ' R (S) " and R2 = Z(S)'/R. 

f 

SUB Rad (CoordO, TempO) 

SHARED NoNodes% 

FOR I - 1 TO NoNodes% 

Rll - Coord (7 , I) * Coord (4, I) 

R12 = Coord (5, I) * Coord(6, I) 

' When computing R2 must allow for fact that R=0 at the first node 
IF I = 1 THEN R2 = Coord(7, I) ELSE R2 = Coord(5, I) / Coord(2 f I) 
Temp (I) = (Rll - R12) A 2 + R2 A 2 
NEXT I 

END SUB 



r ******************************* ROOT ***************************** 
' This subroutine finds the (non-dimensional freq.) A 2 for of the first mode. 
SUB Root (Guess, DT, Cl{), C2(), Answer) 

SHARED NoModes% 


DIM CS (3, 3) 

SELECT CASE NoModes% 

CASE 1 9 Only one term; solve it directly 

Answer = Cl(l, 1) / C2(l, 1) 

IF Answer < 0 THEN PRINT "Eigenvalue is negative!": END 
CASE 2 ' Two terms; use quadratic formula 

RQuad - C2 ( 1 , 1) * C2 (2, 2) - C2(l, 2) * C2(2, 1) 

RMidl = Cl(l, 1) * C2 (2 , 2) + Cl (2, 2) * C2 (1, 1) 

RMid2 = Cl(l, 2) * C2 (2 , 1) + Cl (2, 1) * C2 (1, 2) 

RMid = RMidl - RMid2 

RConst * Cl(l, 1) * Cl (2, 2) - Cl(l, 2) * Cl(2, 1) 

RootTerm *= RMid A 2 - 4 * RQuad * RConst 

IF RootTerm < 0 THEN PRINT "Eigenvalue is complex!": END 
Answerl = (-RMid + SQR (RootTerm) ) / (2 * RQuad) 

Answer2 = (-RMid - SQR (RootTerm) ) / (2 * RQuad) 

IF Answerl < 0 AND Answer2 < 0 THEN PRINT "Eigenvalue is < 0": END 

IF Answerl < 0 THEN Answer = Answer2 

IF Answer2 < 0 THEN Answer - Answerl 

IF Answerl > 0 AND Answer2 > 0 THEN 

IF Answerl < Answer2 THEN Answer = Answerl ELSE Answer = Answer2 
END IF 

CASE 3 ' Three terms; expand the determinant and iterate 

WHILE DT > .000001 

FOR I = 1 TO 3 ' Calculate determinant for "Guess" 

FOR J = 1 TO 3 


CS(I, J) 

= 

C1(I, J) 

+ Guess * C2(I, 

J) 




NEXT J 
NEXT I 

Terml = CS(1, 

i) 

★ 

(CS (2, 

2) * CS (3, 3) - 

CS (3, 

2) * 

CS (2, 

3) ) 

Term2 = CS (1, 

2) 

★ 

(CS (2, 

1) * CS(3, 3) - 

CS (3, 

1) * 

CS (2, 

3) ) 

Term3 = CS (1, 

3) 

★ 

(CS (2 , 

1) * CS (3, 2) - 

CS (3, 

1) * 

CS (2, 

2) ) 

StartVal = Terml 

- 

Term2 

+ Term3 ' 

Determinant 

value 



TestVal = StartVal 

Test = StartVal * TestVal: K = 1 

WHILE Test > 0 9 When Test < 0 we have bracketed the root 

Guess = Guess + DT 
FOR I = 1 TO 3 
FOR J = 1 TO 3 

CS(I f J) = Cl (I, J) + Guess * C2 { I, J) 

NEXT J 


NEXT I 


Terml = CS (1, 

i) 

* (CS (2, 

2) * CS (3, 

3) 

- CS (3, 

2) 

★ 

CS (2, 

3) ) 

Term2 = CS (1, 

2) 

* (CS (2, 

1) * CS (3, 

3) 

- CS (3, 

1) 

★ 

CS (2, 

3)) 

Term3 = CS (1, 

3) 

* (CS (2, 

1) * CS (3, 

2) 

- CS (3, 

1) 

★ 

CS (2, 

2) ) 

TestVal = Terml 

- Term2 + 

Term3 







Test = StartVal 

* TestVal 


K = 

K + 1 





IF K > 100 THEN 

PRINT ”No 

Convergence" : 

END 





WEND 










Guess = Guess - 

DT 

' back 

up one value 

to 

get new 

StartValue 


DT = DT / 10 


' decrease the iteration jump 






WEND 

Answer = Guess 


END SELECT 
END SUB 



9 ****************************** SLOSHM **************************** 

' This subroutine computes the non-dimensional slosh mass of the 
' pendulum model and makes it non-dimensional by dividing by the liquid mass 
' that will completely fill the tank = pi * R0 A 3*LiqVolume/ (FillLevel/100) . 

' The subroutine also computes the ratio of the moment to the force. 

' The mass is returned as Ansi, the ratio is returned as Ans2. 

/ 

SUB SloshM (SO, W(), PS(), PW(), TF ( ) , BI(), SlsFrq, Ansi, Ans2) 

SHARED NoNodes%, NoModes%, FillLevel, BondNo, LiqVolume 
SHARED Ang, ConstMass, ExpMass, LastDZsurf, DeltaElev, SurfElev 

' SO = coordinate matrix for free surface 

9 W() = coordinate matrix for wall 

9 PS ( ) = potential values at free surface nodes 

' PW() = potential values at wall nodes 

' TF ( ) = matrix of trial (eignevalues ) A 2 

' BI() = matrix of T, b” modal participation factors 

' First, compute the nondimensional force and moment amplitudes 

9 The first component is the surface tension effect at the tank walls 

f Sine of the wall angle X at the contact line: 

NN% = NoNodes%: NN1% = NN% - 1 

DeltaZ = W (3, NN% ) - W(3, NN1%) 

DeltaR = W(2, NN% ) - W(2, NN1%) 

SineX = DeltaZ / SQR (DeltaZ A 2 + DeltaR A 2) 

9 Curvature of free surface at the contact line: 

WallCurvl - S (7 , NN%) * S(4, NN%) - S(5, NN%) * S(6, NN%) 

WallCurv2 = S(5, NN% ) / S(2, NN%) 

WallCurv = WallCurvl + WallCurv2 

' Value of the normal derivative of the potential at the contact line: 

WallMass = ConstMass - 1: ' this is parameter "e M 

PotDeriv - 0 

FOR I * 1 TO NoModes% 

PotDeriv - PotDeriv + PS ((I + 1), NN%) * TF(I) * BI(I) 

NEXT I 

PotDeriv = PotDeriv * WallMass * (1 + BondNo) 

' Value of wave amplitude at the contact line: 

WaveAmp - PotDeriv / SQR (SlsFrq) 

' Average the finite element contact angle and the input contact angle: 

RS = S (2, NN%) : RSI - S(2, (NN% - 1)) 

Rw - W (2 , NN%) : RW1 = W(2, (NN% - 1)) 

ZW * W (3, NN% ) : ZW1 = W (3, (NN% - 1) ) 

Terml = (RS - RSI) * (Rw - RW1) + (LastDZsurf) * (ZW - ZW1) 

Term2 = SQR((RS - RSI) A 2 + (LastDZsurf) A 2) 

Term3 = SQR( (Rw - RW1) A 2 + (ZW - ZW1) A 2) 

Cosine * Terml / (Term2 * Term3) 

Sine - SQR (1 - Cosine A 2) 

Angl = ATN(Sine / Cosine) 

AngleTrue = (Angl + Ang * 3.14159 / 180) / 2 
Cosine = COS (AngleTrue) 

Sine = SIN (AngleTrue) 

' The surface tension force: 

Rc = S (2, NN%) : ScriptRc = 1 / WallCurv 

SurfForce * (Cosine - (Rc / ScriptRc) * SineX) * 

SurfForce - SurfForce * WaveAmp 
IF SurfForce > 0 THEN SurfForce = 0 


Cosine / Sine 



' The surface tension moment: 

Fc = (W(3, NoNodes%) - W(3, 1)) - SurfElev 

FirstPart = Cosine - ( (Rc / ScriptRc) + (Rc / Fc) * Sine) * SineX 
FirstPart = FirstPart * Cosine / Sine 
SurfMoment = FirstPart * Fc * WaveAmp 

' Pressure component of force and moment: 

PresForce = 0: PresMoment = 0 

FOR I = 1 TO NoNodes% 

SELECT CASE I ' get the right DeltaZ 

CASE 1 

DZ = (W (3, 2) - W{3, 1) ) / 2 
CASE NoNodes% 

DZ « (W (3, NN%) - W{3, (NN% -1))) / 2 
CASE ELSE 

DZ = (W(3, (I + 1)) - W (3, (I - 1))) / 2 
END SELECT 

FOR J = 1 TO NoModes% 

Rw = W(2, I) 

Z = DeltaElev + W(3, I) — SurfElev ' referred to free surface 

TrialPotential = PW((J +1), I) * BI ( J) 

PresForce = PresForce + TrialPotential * Rw * DZ 
PresMoment = PresMoment + TrialPotential * Rw * z * DZ 
NEXT J 
NEXT I 


Force = PresForce + (SurfForce) / ((1 + BondNo) * SQR(SlsFrq)) 
Moment = PresMoment + (SurfMoment) / ( (1 + BondNo) * SQR(SlsFrq)) 

IF Force < 0 THEN CLS : PRINT "The Slosh Force is Negative!": END 
Ans2 = Moment / Force 


' The non-dimensional kinetic energy of the liquid: 

KinEnergy = 0 
SMax = S (3, NN% ) 

FOR I = 1 TO NoNodes% 

SELECT CASE I 

CASE 1 ' get the right value ofDeltaS 

DS = (S (3, 2) - S(3, 1) ) / 2 
CASE NoNodes% 


DS = (S (3, NN%) - S (3, (NN% - 1))) / 2 
CASE ELSE 

DS = (S (3, (I + 1)) - S (3, (I - 1))) / 2 
END SELECT 

RS = S (2, I) : SS = S (3, I) 

SurfMass = (ConstMass - (SS / SMax) * ExpMass) * (1 + BondNo) 
PotDeriv = 0 


Pot = 0 

FOR J - 1 TO NoModes% ' dPhi/dN at the integration point 

PotDeriv = PotDeriv + BI ( J) * TF(J) * SurfMass * PS ( ( J + 1) , I) 
NEXT J 

FOR J = 1 TO NoModes% ' Phi at the integration point 

Pot = Pot + PS((J + 1), I) * BI (J) 

NEXT J 

KinEnergy = KinEnergy + PotDeriv * Pot * RS * DS 

NEXT I 


' Now compute the slosh mass. The mass is ratioed to the liquid mass that 
' fills the tank. 

Ansi = (Force) A 2 * (FillLevel / 100) / (KinEnergy * LiqVolume) 

END SUB 



r a**************************** TEXT IN ********************************** 

' This subroutine is taken from the QUICKPAK set of functions. It is a 
' data entry routine that allows the cursor to be moved all over the 
' screen and permits the entered data to be corrected at any time. 

' Entry Parameters 

' t$ = string to be input or edited (use the name of your choice) 

' Max = maximum number of characters allowed 

' Coir is the combined foreground and background colors that will be used 
' CapsOn = force automatic conversion to upper case if 1 
' Note: CapsOn is not used here 

' NumOnly = allow only numeric input if 1 

' Exit Parameters 

' T$ = final edited string (the name passed as input to the function) 

' ExitCode indicates how editing was terminated - 

' 0 = Enter, Tab, Down-Arrow, Right-Arrow past end, or field filled 

' 1 = Shift-Tab, Up-Arrow, or Left-Arrow past beginning 

' 2 = Esc key pressed 

' Local Variables 

' x$ is a copy of the string while being input or edited 

' Insert holds status of insert mode 

' Curpo holds current cursor position relative to the beginning of the line 
' Length keeps track of the current length of the string 
' cir = 1 if the monitor is a color monitor, for setting cursor size 
' A and A$ are temporary scratch variables 

SUB Textlnl (T$ , Max, NumOnly, CapsOn, ExitCode, Coir) STATIC 
DEFINT A-Z 
TInitialize : 

Clr = 0 

IF Peekl% (0, &H463) <> &HB4 THEN Clr = 1 
X$ = T$ 

TC: 

ExitCode = 0: Insert = 0: Curpo * 1 

Length = LEN(X$) 

IF Length > Max THEN EXIT SUB 

X$ = X$ + SPACE$ (Max - Length) 

QPrint X$, Coir, -1 
LOCATE , , 1 

GOSUB TInsertOff 

TGetKey : 

' disallow insert if cursor past end 
IF Curpo > Length AND Insert <> 0 THEN GOSUB TInsertOff 

IF Curpo > Max GOTO TEnter 'field is filled, handle as Enter key 

A$ = INKEY $ 

IF A$ = GOTO TGetKey 
IF LEN (A$ ) - 1 GOTO TRegularKey 

A$ = RIGHT$(A$, 1) 


'already to big to edit 
'pad with trailing spaces 

'turn on the cursor 
'set cursor size according to display 


'determine monitor type 

' X$ is a working copy of 
' the input string 

'initialize flags 


'it was an extended key, get the code 



ON INSTR(CHR$ (15) + " . GHKMOPRSstu" + CHR$(19), A$) GOTO TShiftTab, TClear 
THome, TUp, TLeft, TRight, TEndKey, TDown, Tins, TDel, TCtrlLeft, TCtrlRight 
TCtrlEnd, TRestore 


GOTO TGetKey 

TShiftTab: 

ExitCode = 1 
GOTO TEnter 

TClear : 

X$ = "" 

GOSUB TInsertOff 

LOCATE , POS(O) - (Curpo - 1) 

GOTO TC 

THome : 

LOCATE , POS(O) - (Curpo - 1) 
Curpo = 1 
GOTO TGetKey 

TUp: 

ExitCode = 1 
GOTO TEnter 

TLeft : 

IF Curpo = 1 GOTO TShiftTab 

Curpo = Curpo - 1 
LOCATE , POS ( 0 ) - 1 
GOTO TGetKey 

TRight : 

Curpo = Curpo + 1 
LOCATE , POS(O) + 1 
GOTO TGetKey 

TEndKey : 

LOCATE , POS(O) + (Length - Curpo) 
Curpo = Length + 1 
GOTO TGetKey 

TDown : 

GOTO TEnter 

Tins: 


'none of the above, get again 


'user wants to go back a field 
'handle as if it were the Enter key 

'Alt-C, erase the current string 
'clear insert mode and restore cursor 

'and start all over again 


'put cursor at beginning of line 
'show cursor as being on 1st character 


'user wants to go back a field 
'handle as if it were the Enter key 


'cursor is on the first character, 
'handle as if it were a Shift-Tab 
'update Curpo 
'back up the cursor 


'update Curpo 

'advance the cursor on the screen 


+ 1 'put cursor at the end of the line 
'update Curpo 


IF Insert THEN 'insert is already on, turn it off 

GOSUB TInsertOff 
GOTO TGetKey 
END IF 


IF Curpo > Length GOTO TGetKey 
IF Length « Max GOTO TGetKey 

Insert = 1 
IF Clr THEN 

LOCATE , , , 0, 7 
ELSE 

LOCATE , , , 0, 13 
END IF 

GOTO TGetKey 
TDel : 


' ignore Ins if cursor is past the end 
'also ignore if field is full 

'set the insert flag 

'set cursor size according to display 


IF Curpo > Length GOTO TGetKey 


ignore Del if cursor is past end 



'slide all characters left one position, 
MID$(X$, Curpo) = MID$(X$, Curpo + 1) 
QPrint MID$(X$, Curpo), Coir, -1 

Length = Length - 1 
GOTO TGetKey 

TCtrlLef t : 

IF Curpo = 1 GOTO TGetKey 
A = Curpo 

'we're within a word, find beginning 
IF MID$(X$, Curpo - 1, 1) <> " 

TSeekLeftl : 

IF Curpo = 1 GOTO TCtrlLeftExit 

IF MID$(X$, Curpo - 1, 1) = " " 

Curpo = Curpo - 1 
GOTO TSeekLeftl 
END IF 

TSeekLef t2 : 

IF Curpo = 1 GOTO TCtrlLeftExit 
IF MID$(X$, Curpo - 1, 1) <> " 

Curpo = Curpo - 1 
GOTO TSeekLeft2 
END IF 

TCtrlLeftExit : 

LOCATE , POS(O) - (A - Curpo) 

GOTO TGetKey 

TCtrlRight : 

A = Curpo 

TSeekRightl : 

IF A > Length GOTO TGetKey 

IF MID$(X$, A, 1) <> " " THEN 
A = A + 1 
GOTO TSeekRightl 
END IF 

TSeekRight2 : 

IF A > Length GOTO TGetKey 

IF MID$(X$, A, 1) = " " THEN 
A = A + 1 
GOTO TSeekRight2 
END IF 

LOCATE , POS(O) + (A - Curpo) 

Curpo = A 
GOTO TGetKey 

TCtrlEnd: 

IF Curpo > Length GOTO TGetKey 

QPrint SPACES (Length - Curpo + 1), 
Length = Curpo - 1 
GOTO TGetKey 


add a trailing space and re-print 

M N 

'show string as one character shorter 

'at the beginning, ignore 
'save cursor position 

GOTO TSeekLeft2 


'at the beginning, give up 

THEN 

' seek previous non-blank character 


'at the beginning, give up 

THEN 

'seek character preceeded by a blank 


'position the cursor 


'save cursor position 

'at the end, give up 

'consider next character 
' seek next blank space 


'at the end, give up 

'consider next character 
'seek next non-blank character 

'position the cursor 

'show cursor as being on the next word 
'get another keypress 

'cursor is past the end, ignore 

Coir, -1 'blank from cursor to end 
' show the length being at the cursor 
'get another keypress 


TRestore : 



LOCATE , POS(O) - (Curpo - 1) 
GOTO TInitialize 


'locate cursor at beginning of line, 

' and start all over again 

TRegularKey : 

IF A$ < " " THEN ' a control key 

ON INSTR (CHR$ ( 8 ) + CHR$<9) + CHR$(13) + CHR$<27), A$) GOTO TBackspace, 
TTabKey, TEnter, TEscape 

GOTO TGetKey 'none of the above 

END IF 


IF CapsOn THEN 'convert to upper case if requested 

IF A$ >= "a" AND A$ <= "z" THEN A$ = CHR$(ASC(A$) AND 95) 

END IF 


IF NumOnly THEN 'disallow non-numeric if requested 

IF A$ < "0" OR A$ > "9" THEN 
PLAY "LI 603EC” 

GOTO TGetKey 
END IF 
END IF 


QPrint A$, Coir, -1 'print character 

LOCATE f POS(O) + 1 

Curpo = Curpo + 1 ' show cursor being ahead 

IF Insert GOTO THandlelnsert 

MID$ (X$, Curpo - 1, 1) = A$ 'assign the character 

'cursor is past end, increase length 

IF Curpo > Length + 1 THEN Length = Curpo - 1 

'field complete, handle as Enter key 

IF Length = Max AND Curpo > Length GOTO TEnter 

GOTO TGetKey 

THandlelnsert : 


Length = Length + 1 'show string being 1 character longer 

MID$(X$, Curpo) = MID$(X$, Curpo - 1) 'move characters one position ahead 

MID$(X$, Curpo - 1, 1) = A$ 'assign the current character 

QPrint MID$(X$, Curpo, Length - Curpo + 1), Coir, -1 're-print X$ 

IF Length — Max GOTO TEnter 'field complete, handle as Enter key 

GOTO TGetKey 

TBackspace : 


IF Curpo * 1 GOTO TGetKey 
Curpo * Curpo - 1 
LOCATE , POS(O) - 1 
GOTO TDel 

TTabKey : 

TEnter : 


'can't back up any more, ignore 
' show cursor being 1 character before 
'back up the cursor 

'handle as if it were the Delete key 

'reserved for your Tab routine if you 
' want to handle it differently 


GOSUB TInsertOff 

X$ = LEFT$(X$, Length) 

T$ = X$ 

LOCATE , , 0 
EXIT SUB 


'clear insert, restore cursor size 
'retain only the current length 

'assign the string 
'turn off the cursor 


TEscape : 

ExitCode = 2 
GOTO TEnter 


TInsertOff : 


' show that the user pressed Escape 
'handle as if it were the Enter Key 



'clear Insert mode and restore cursor, depending on monitor type 

Insert * 0 
IF Clr THEN 

LOCATE , , , 6, 1 

ELSE 

LOCATE , , , 12, 13 
END IF 
RETURN 


END SUB 



a******************************* VECTOR ************************** 

This subroutine computes the eigenvectors or modal particpation factors 


SUB Vector (F2, CIO, C2(), Result {) ) 

DEFSNG A-Z 

SHARED NoModes% 

DIM CR (3, 3) 

SELECT CASE NoModes% 

CASE 1 

Result (1) = 1 
CASE 2 

Result (1) = 1 

Result (2) = -(Cl(l, 1) + F2 * C2(l, 1)) / (Cl (1, 2) + F2 * C2(l, 2)) 
CASE 3 

Result (1) = 1 
FOR I = 1 TO 3 
FOR J = 1 TO 3 

CR ( I , J) = C1(I, J) + F2 * C2 ( I, J) 

NEXT J 
NEXT I 

TermTop = CR(2, 2) * CR(1, 1) - CR(2, 1) * CR(1, 2) 

TermBot = CR(1, 2) * CR(2, 3) - CR(1, 3) + CR(2, 2) 

Result (3) = TermTop / TermBot 

Result (2) = (-CR ( 1, 3) * Result (3) - CR(1, 1)) / CR(1, 2) 

END SELECT 

END SUB 



r *★***★*★★**★*★*★*★*★**★**★*** WAVE HITE **★* ********************* * 

' This subroutine computes the wave height of the trial functions. 

' However, the frequency multiplication term is not included. 

r 

SUB WaveHite (Coord (), ModeS 0, TempNoO) 

SHARED NoModes%, NoNodes%, ConstMass, ExpMass 

SMax = Coord (3, NoNodes%) ' Value of S coord at contact point 

FOR I - 1 TO NoNodes% 

DistMass - ConstMass - (Coord (3, I) / SMax) A ExpMass 
FOR J » 1 TO NoModes% 

TempNo(J, I) * DistMass * ModeS ( (J +1), I) 

NEXT J 
NEXT I 

END SUB 



